function y_non = fiscal_reaction(x_,T_,S_)

%% Choose Fiscal Reaction

% S_.options.fiscal = 'gss_tauK'; 

% 'gss_tr': government maintains fixed level of government spending and lump-sum transfers adjust to satisfy budget (this is the baseline)
% 'gss_b': government maintains fixed level of government spending and public debt adjusts to satisfy budget
% 'gss_tauc': government maintains fixed level of government spending and consumption tax rate adjusts to satisfy budget
% 'gss_tauw': government maintains fixed level of government spending and labour income tax rate adjusts to satisfy budget
% 'gss_tauK': government maintains fixed level of government spending and capital income tax rate adjusts to satisfy budget
% 'gss_taxprop': government maintains fixed level of government spending and all tax rates adjusts proportionally to satisfy budget

% 'gyss_tr': government maintains fixed government spending-to-GDP ratio and lump-sum transfers adjust to satisfy budget
% 'gyss_b': government maintains fixed government spending-to-GDP ratio and public debt adjusts to satisfy budget
% 'gyss_tauc': government maintains fixed government spending-to-GDP ratio and consumption tax rate adjusts to satisfy budget
% 'gyss_tauw': government maintains fixed government spending-to-GDP ratio and labour income tax rate adjusts to satisfy budget
% 'gyss_tauK': government maintains fixed government spending-to-GDP ratio and capital income tax rate adjusts to satisfy budget
% 'gyss_taxprop': government maintains fixed government spending-to-GDP ratio and all tax rates adjusts proportionally to satisfy budget

% 'gyss_tyss_b': government maintains fixed government spending-to-GDP ratio and lump-sum transfert-to-GDP and public debt adjusts to satisfy budget
% 'gyss_tyss_tauc': government maintains fixed government spending-to-GDP ratio and lump-sum transfert-to-GD and consumption tax rate adjusts to satisfy budget
% 'gyss_tyss_tauw': government maintains fixed government spending-to-GDP ratio and lump-sum transfert-to-GD and labour income tax rate adjusts to satisfy budget
% 'gyss_tyss_tauK': government maintains fixed government spending-to-GDP ratio and lump-sum transfert-to-GD and capital income tax rate adjusts to satisfy budget
% 'gyss_tyss_taxprop': government maintains fixed government spending-to-GDP ratio and lump-sum transfert-to-GD and all tax rates adjusts proportionally to satisfy budget


%% Data and parameters

aus_data;

Tz = T_(1);
Tsig = T_(2);

%% Parameters

S_.param = param_func(x_);

%% Variable and dimensions
    
v_ = var_dim ;
    
n     = v_.dim.n    ;
n_    = v_.dim.n_    ;
ns    = v_.dim.ns    ;
           

%% Structure 1 (z_ss low, gy_ss)

S_.options.initial = 1;

S_.initvals = [ 0.3895      % ln c -- 1
                -0.6735     % ln invest -- 2
                -0.4527     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3254      % ln k -- 5
                -0.0124     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9550      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3970      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                0.0873      % tr (can be negative) -- 17
                3.1845];    % p_b -- 18         

st_1 = fiscalreaction_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; 


% Set initial steady state fiscal values
S_.param.b_ss = log(st_1.y_steady(15));
S_.param.g_ss = st_1.y_steady(16);
S_.param.tr_ss = log(st_1.y_steady(17));
S_.param.ty_ss = log(st_1.y_steady(17))./st_1.y_steady(9);

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

S_.options.initial = 0;
S_.param.z_ss = S_.param.z_ss + S_.param.Delta_z ; 

%%%%%%%%%%%%%%%% Case baseline %%%%%%%%%%%%%%%%
if strcmp('baseline',S_.options.fiscal)

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_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -1.4923     % ln gy -- 16
                0.0568  % tr (can be negative) -- 17
                3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case baseline %%%%%%%%%%%%%%%%
elseif strcmp('gyss_baseline',S_.options.fiscal)

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_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                0.0568    % tr (can be negative) -- 17 
                3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case Short Term Debt %%%%%%%%%%%%%%%%
elseif strcmp('short_term_debt',S_.options.fiscal)

S_.param = param_func(x_);
S_.param.kappa_b = 0;

S_.options.initial = 1;

S_.initvals = [ 0.3895      % ln c -- 1
                -0.6735     % ln invest -- 2
                -0.4527     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3254      % ln k -- 5
                -0.0124     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9550      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3970      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                0.0873      % tr (can be negative) -- 17
                3.1845];    % p_b -- 18         

st_1 = fiscalreaction_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; 

S_.options.initial = 0;

S_.param.z_ss = S_.param.z_ss + S_.param.Delta_z ; 

% Set initial steady state fiscal values
S_.param.b_ss = log(st_1.y_steady(15));
S_.param.g_ss = st_1.y_steady(16);
S_.param.tr_ss = log(st_1.y_steady(17));
S_.param.ty_ss = log(st_1.y_steady(17))./st_1.y_steady(9);

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_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -1.4923     % ln gy -- 16
                0.0568  % tr (can be negative) -- 17
                3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);


%%%%%%%%%%%%%%%% Case gss_tr %%%%%%%%%%%%%%%%
elseif strcmp('gss_tr',S_.options.fiscal)

S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -1.4923     % ln gy -- 16
                0.0568  % tr (can be negative) -- 17
                3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gss_b %%%%%%%%%%%%%%%%
elseif strcmp('gss_b',S_.options.fiscal)

S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                6.7890      % b  (can be negative) -- 15
                -1.4923     % ln gy -- 16
                2.6288      % by (can be negative) -- 17
                3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gss_tauc %%%%%%%%%%%%%%%%
elseif strcmp('gss_tauc',S_.options.fiscal)

S_.initvals = [ 0.4100      % ln c -- 1
                -0.6741     % ln invest -- 2
                -0.4539     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3248      % ln k -- 5
                -0.0130     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9544      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3961      % b  (can be negative) -- 15
                -1.4923     % ln gy -- 16
                -3.0724     % tau_c
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gss_tauw %%%%%%%%%%%%%%%%
elseif strcmp('gss_tauw',S_.options.fiscal)

S_.initvals = [ 0.4104      % ln c -- 1
                -0.6738     % ln invest -- 2
                -0.4734     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3250      % ln k -- 5
                -0.0127     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9547      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3961      % b  (can be negative) -- 15
                -1.4982     % ln gy -- 16
                -1.8190    % tau_w -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gss_tauK %%%%%%%%%%%%%%%%
elseif strcmp('gss_tauK',S_.options.fiscal)

S_.initvals = [ 0.4199      % ln c -- 1
                -0.5918     % ln invest -- 2
                -0.4828     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.4071      % ln k -- 5
                -0.0151     % ln l -- 6
                0.6494      % ln w -- 7
                -3.6663     % ln r_k -- 8
                0.9767      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.4276      % b  (can be negative) -- 15
                -1.5203     % ln gy -- 16
                -1.3609   % tau_K -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gss_taxprop %%%%%%%%%%%%%%%%
elseif strcmp('gss_taxprop',S_.options.fiscal)

S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191     % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -1.4923     % ln gy -- 16
                -2.9724   % prop -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tr %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tr',S_.options.fiscal)

S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                0.0568    % tr (can be negative) -- 17 
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_b %%%%%%%%%%%%%%%%
elseif strcmp('gyss_b',S_.options.fiscal)
    
S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                6.7890      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                2.6288    % by (can be negative) -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tauc %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tauc',S_.options.fiscal)

S_.initvals = [ 0.4100      % ln c -- 1
                -0.6741     % ln invest -- 2
                -0.4539     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3248      % ln k -- 5
                -0.0130     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9544      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3961      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                -3.0724    % tau_c -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tauw %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tauw',S_.options.fiscal)

S_.initvals = [ 0.4104      % ln c -- 1
                -0.6738     % ln invest -- 2
                -0.4734     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3250      % ln k -- 5
                -0.0127     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9547      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3961      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                -1.8190  % tau_w -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tauK %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tauK',S_.options.fiscal)

S_.initvals = [ 0.4199      % ln c -- 1
                -0.5918     % ln invest -- 2
                -0.4828     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.4071      % ln k -- 5
                -0.0151     % ln l -- 6
                0.6494      % ln w -- 7
                -3.6663     % ln r_k -- 8
                0.9767      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.4276      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                -2.4609   % tau_K -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_taxprop %%%%%%%%%%%%%%%%
elseif strcmp('gyss_taxprop',S_.options.fiscal)

S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191     % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -0.4923     % ln gy -- 16
                -5.0672     % prop -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tyss_b %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tyss_b',S_.options.fiscal)
    
S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191      % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                6.7890      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                2.6288    % by (can be negative) -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tyss_tauc %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tyss_tauc',S_.options.fiscal)

S_.initvals = [ 0.4100      % ln c -- 1
                -0.6741     % ln invest -- 2
                -0.4539     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3248      % ln k -- 5
                -0.0130     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9544      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3961      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                -3.0724   % tau_c -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tyss_tauw %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tyss_tauw',S_.options.fiscal)

S_.initvals = [ 0.4104      % ln c -- 1
                -0.6738     % ln invest -- 2
                -0.4734     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3250      % ln k -- 5
                -0.0127     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9547      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3961      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                -1.8190   % tau_w -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tyss_tauK %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tyss_tauK',S_.options.fiscal)

S_.initvals = [ 0.4199      % ln c -- 1
                -0.5918     % ln invest -- 2
                -0.4828     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.4071      % ln k -- 5
                -0.0151     % ln l -- 6
                0.6494      % ln w -- 7
                -3.6663     % ln r_k -- 8
                0.9767      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.4276      % b  (can be negative) -- 15
                -0.4900     % ln g -- 16
                -2.4609   % tau_K -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);

%%%%%%%%%%%%%%%% Case gyss_tyss_taxprop %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tyss_taxprop',S_.options.fiscal)

S_.initvals = [ 0.4022      % ln c -- 1
                -0.6798     % ln invest -- 2
                -0.4653     % ln lambda -- 3
                0.0069      % ln r_F -- 4
                3.3191     % ln k -- 5
                -0.0186     % ln l -- 6
                0.6249      % ln w -- 7
                -3.6063     % ln r_k -- 8
                0.9488      % ln y --9
                0.0000      % nx (can be negative) -- 10
                0.0000      % ca (can be negative) -- 11
                0.0020      % r_star -- 12
                0.0000      % b_F (can be negative) -- 13
                0.0069      % ln r -- 14
                1.3882      % b  (can be negative) -- 15
                -0.4923     % ln gy -- 16
                -5.0672   % prop -- 17  
                 3.1845];    % p_b -- 18

st_2 = fiscalreaction_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_,:);


end


%% Plot simulated data against data

T = size(Zobs,2);

 T_extend = 200; 

 T = size(Zobs,2) + T_extend;

y_simul = zeros(n,T); 
y_simul_bp = zeros(n_,T);

y_simul_break = zeros(n,T); 
y_non = zeros(n,T); 

OMEGA = st_1.OMEGA;

% Set Random number generator (so that the results are the same each time)
rng('default')
rng(171);
RandStream.getGlobalStream;

e_(1:ns,1:Tsig-1) = (S_.param.mu_sig*OMEGA).^.5 * randn(ns,Tsig-1);
e_(1:ns,Tsig:T) = OMEGA.^.5 * randn(ns,T-Tsig+1);

y_simul(:,1) = y_bar_sims_1 + S2_1*e_(:,1);
y_simul_bp(:,1) = y_bar_1 + G_1*e_(:,1);
y_simul_break(:,1) = y_bar_sims_1 + S2_1*e_(:,1);
y_non(:,1) = y_bar_sims_1;

for t = 2:T
    y_simul(:,t) = S0_1 + S1_1 * y_simul(:,t-1) + S2_1*e_(:,t);
    y_simul_bp(:,t) = C_1 + Q_1 * y_simul_bp(:,t-1) + G_1*e_(:,t);
end

% with the break in z_ss and g_ss
for t = 2:Tz-1
    y_simul_break(:,t) = S0_1 + S1_1 * y_simul_break(:,t-1) + S2_1*e_(:,t);
    y_non(:,t) = S0_1 + S1_1*y_non(:,t-1);
end
for t = Tz:T
    y_simul_break(:,t) = S0_2 + S1_2 * y_simul_break(:,t-1) + S2_2*e_(:,t);
    y_non(:,t) = S0_2 + S1_2*y_non(:,t-1);
end

% Order in Zobs
% 1. dy_obs
% 2. dc_obs_adj
% 3. nx_obs
% 4. dw_obs
% 5. r_real_obs
% 6. r_real_star_obs
% 7. gy_obs
% 8. by_obs
% 9. rev_c_y_obs_adj 
% 10. rev_w_y_obs
% 11. rev_K_y_obs

% figure(1)
% subplot(4,3,1),plot(t_,Zobs(1,:),'k',t_,y_simul(v_.dy_obs,:),'r'), title('dy')
% subplot(4,3,2),plot(t_,Zobs(2,:),'k',t_,y_simul(v_.dc_obs,:),'r'), title('dc')
% subplot(4,3,3),plot(t_,Zobs(3,:),'k',t_,y_simul(v_.nx_obs,:),'r'), title('nx')
% subplot(4,3,4),plot(t_,Zobs(4,:),'k',t_,y_simul(v_.dw_obs,:),'r'), title('dw')
% subplot(4,3,5), plot(t_,Zobs(5,:),'k',t_,y_simul(v_.r_obs,:),'r'), title('r')
% subplot(4,3,6), plot(t_,Zobs(6,:),'k',t_,y_simul(v_.r_star_obs,:),'r'), title('r')
% subplot(4,3,7), plot(t_,exp(Zobs(7,:)),'k',t_,exp(y_simul(v_.gy_obs,:)),'r'), title('G/Y')
% subplot(4,3,8), plot(t_,Zobs(8,:),'k',t_,y_simul(v_.by_obs,:),'r'), title('B/Y')
% subplot(4,3,9), plot(t_,Zobs(9,:),'k',t_,y_simul(v_.rev_c_y_obs,:),'r'), title('rev_c_y')
% subplot(4,3,10), plot(t_,Zobs(10,:),'k',t_,y_simul(v_.rev_w_y_obs,:),'r'), title('rev_w_y')
% subplot(4,3,11), plot(t_,Zobs(11,:),'k',t_,y_simul(v_.rev_K_y_obs,:),'r'), title('rev_K_y')
% 
% figure(2)
% subplot(4,3,1),plot(t_,Zobs(1,:),'k',t_,y_simul_break(v_.dy_obs,:),'r'), title('dy')
% subplot(4,3,2),plot(t_,Zobs(2,:),'k',t_,y_simul_break(v_.dc_obs,:),'r'), title('dc')
% subplot(4,3,3),plot(t_,Zobs(3,:),'k',t_,y_simul_break(v_.nx_obs,:),'r'), title('nx')
% subplot(4,3,4),plot(t_,Zobs(4,:),'k',t_,y_simul_break(v_.dw_obs,:),'r'), title('dw')
% subplot(4,3,5), plot(t_,Zobs(5,:),'k',t_,y_simul_break(v_.r_obs,:),'r'), title('r')
% subplot(4,3,6), plot(t_,Zobs(6,:),'k',t_,y_simul_break(v_.r_star_obs,:),'r'), title('r')
% subplot(4,3,7), plot(t_,exp(Zobs(7,:)),'k',t_,exp(y_simul_break(v_.gy_obs,:)),'r'), title('G/Y')
% subplot(4,3,8), plot(t_,Zobs(8,:),'k',t_,y_simul_break(v_.by_obs,:),'r'), title('B/Y')
% subplot(4,3,9), plot(t_,Zobs(9,:),'k',t_,y_simul_break(v_.rev_c_y_obs,:),'r'), title('rev_c_y')
% subplot(4,3,10), plot(t_,Zobs(10,:),'k',t_,y_simul_break(v_.rev_w_y_obs,:),'r'), title('rev_w_y')
% subplot(4,3,11), plot(t_,Zobs(11,:),'k',t_,y_simul_break(v_.rev_K_y_obs,:),'r'), title('rev_K_y')


% figure(3)
% subplot(4,3,1),plot(t_,Zobs(1,:),'k',t_,y_non(v_.dy_obs,:),'r'), title('dy')
% subplot(4,3,2),plot(t_,Zobs(2,:),'k',t_,y_non(v_.dc_obs,:),'r'), title('dc')
% subplot(4,3,3),plot(t_,Zobs(3,:),'k',t_,y_non(v_.nx_obs,:),'r'), title('nx')
% subplot(4,3,4),plot(t_,Zobs(4,:),'k',t_,y_non(v_.dw_obs,:),'r'), title('dw')
% subplot(4,3,5), plot(t_,Zobs(5,:),'k',t_,y_non(v_.r_obs,:),'r'), title('r')
% subplot(4,3,6), plot(t_,Zobs(6,:),'k',t_,y_non(v_.r_star_obs,:),'r'), title('r')
% subplot(4,3,7), plot(t_,exp(Zobs(7,:)),'k',t_,exp(y_non(v_.gy_obs,:)),'r'), title('G/Y')
% subplot(4,3,8), plot(t_,Zobs(8,:),'k',t_,y_non(v_.by_obs,:),'r'), title('B/Y')
% subplot(4,3,9), plot(t_,Zobs(9,:),'k',t_,y_non(v_.rev_c_y_obs,:),'r'), title('rev_c_y')
% subplot(4,3,10), plot(t_,Zobs(10,:),'k',t_,y_non(v_.rev_w_y_obs,:),'r'), title('rev_w_y')
% subplot(4,3,11), plot(t_,Zobs(11,:),'k',t_,y_non(v_.rev_K_y_obs,:),'r'), title('rev_K_y')



