function st_mat = fiscalreaction_matrices(S_)

%% This function takes structural parameters in S_ for the small open economy real business
% cycle model and produces the structural matrices in Sims form.


% Firms

alpha = S_.param.alpha;    % Cobb-Douglas production function parameter
z_ss  = S_.param.z_ss;     % growth rate of labour-augmenting TFP
A_ss  = S_.param.A_ss;     % steady state level of stationary TFP

% Households
beta    = S_.param.beta  ;  % household's discount factor
delta   = S_.param.delta ;   % capital depreciation
gamma   = S_.param.gamma ;    % labour supply 
nu      = S_.param.nu ;       % labour supply (Frisch)
Upsilon = S_.param.Upsilon ; % 3.32;    % Upislon'' investment adjustment cost
h       = S_.param.h ; % habits

% Foreign sector
zeta_b_ss = S_.param.zeta_b_ss;   
S_.param.rstar_ss = S_.param.z_ss / S_.param.beta*exp(-S_.param.zeta_b_ss); % steady state of foreign interest rate 
by_F_ss  = S_.param.by_F_ss ;    % net foreign asset to output ratio (+ asset / (-) debt)
psi_b    = S_.param.psi_b ;      % governs country-risk premium to net foreign asset position

% Government
by_ss    = S_.param.by_ss ; % steady state debt to output ratio
gy_ss    = S_.param.gy_ss ;   % steady state government spending to output ratio
%g_ss     = S_.param.g_ss ;  % steady state government spending per unit of effective labour
tau_c_ss = S_.param.tau_c_ss ;  % steady state tax rate on consumption
tau_w_ss = S_.param.tau_w_ss ;   % steady state tax rate on labour income
tau_K_ss = S_.param.tau_K_ss ;   % steady state tax rate on corporate income

%  -- Fiscal reaction function parameters
rho_g  = S_.param.rho_g ;   % persistence of government spending
rho2_g  = S_.param.rho2_g ;   % persistence of government spending
gam_gb = S_.param.gam_gb ; % response of government spending to b/y

rho_c  = S_.param.rho_c ;  % persistence of tax rate on consumption
rho2_c  = S_.param.rho2_c ;  % persistence of tax rate on consumption
gam_cb = S_.param.gam_cb;  % response of tax rate on consumption to b/y

rho_w  = S_.param.rho_w;  % persistence of tax rate on labour income
rho2_w  = S_.param.rho2_w ;  % persistence of tax rate on labour income
gam_wb = S_.param.gam_wb ;  % response of tax rate on labour income to b/y

rho_K  = S_.param.rho_K ;   % persistence of tax rate on corporate income
rho2_K  = S_.param.rho2_K ;  % persistence of tax rate on corporate income
gam_Kb = S_.param.gam_Kb;  % response of tax rate on corporate income to b/y

rho_tau  = S_.param.rho_tau ; % persistence of lump-sump transfers
rho2_tau  = S_.param.rho2_tau ; % persistence of lump-sump transfers
gam_taub = S_.param.gam_taub ; % response of transfers to b/y

%  -- automatic stabilizers
gam_gy    =  S_.param.gam_gy  ; % response of government spending to y
gam_cy    =  S_.param.gam_cy  ; % response of tax rate on consumption to y
gam_wy    =  S_.param.gam_wy  ; % response of tax rate on labour income to y
gam_Ky    =  S_.param.gam_Ky  ; % response of tax rate on corporate income to y
gam_tauy  =  S_.param.gam_tauy; % response of transfers to y

% AR exogenous shock processes

rho_z      = S_.param.rho_z ;      % Persistence of TFP growth
rho_r_star = S_.param.rho_r_star;  % Persistence of foreign interest rate
rho_zeta   = S_.param.rho_zeta;    % Persistence of intertemporal preference shock
rho_L      = S_.param.rho_L ;      % Persistence of labour supply shock
rho_I      = S_.param.rho_I ;      % Persistence of investment shock
rho_b      = S_.param.rho_b;  % Persistence of risk premium shock
rho_vg     = S_.param.rho_vg;      % Persistence of government spending shock
rho_by     = S_.param.rho_by ;     % Persistence of time-varying by ratio

% Standard deviation of shocks
%  - Fiscal policy shocks
sig_g    = S_.param.sig_g ;
sig_tauc = S_.param.sig_tauc;
sig_tauw = S_.param.sig_tauw;
sig_tauK = S_.param.sig_tauK ;
sig_tau  = S_.param.sig_tau ;

%  - Structural shocks
sig_z     = S_.param.sig_z; 
sig_rstar = S_.param.sig_rstar;
sig_zeta  = S_.param.sig_zeta ;
sig_L     = S_.param.sig_L ;
sig_I     = S_.param.sig_I ;
sig_b     = S_.param.sig_b;

sigma = S_.param.sigma;
kappa_b = S_.param.kappa_b;

%% Calculates steady state

max_iter_steady = 100;

if S_.options.initial

    [y_steady, Out_ ] = steadystate_rbc(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tr_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_); 
    phi_ss    = lambda_ss;
    tax_ss    = g_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

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

    [y_steady, Out_ ] = steadystate_rbc_gss(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss      = y_steady(l_);     l_ = l_ + 1;
    tr_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_); 
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss +  ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss      = y_steady(l_);     l_ = l_ + 1;
    tr_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_); 
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss +  ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tr(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tr_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss_tr(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss      = y_steady(l_);     l_ = l_ + 1;
    tr_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_); 
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss +  ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss_b(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss     = y_steady(l_);      l_ = l_ + 1;
    by_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_); 
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;
    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss_tauc(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss     = y_steady(l_);      l_ = l_ + 1;
    tau_c_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_); 
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;
    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss_tauw(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss     = y_steady(l_);      l_ = l_ + 1;
    tau_w_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;
    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss_tauK(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss     = y_steady(l_);      l_ = l_ + 1;
    tau_K_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;
    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gss_taxprop(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    gy_ss     = y_steady(l_);      l_ = l_ + 1;
    prop      = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    g_ss = S_.param.g_ss;
    tr_ss = S_.param.tr_ss;

    tau_c_ss = tau_c_ss*(1-prop);
    tau_w_ss = tau_w_ss*(1-prop);
    tau_K_ss = tau_K_ss*(1-prop);

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tr(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tr_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_b(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    by_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tauc(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tau_c_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tauw(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tau_w_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tauK(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tau_K_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.tr_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_taxprop(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    prop      = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.tr_ss;

    tau_c_ss = tau_c_ss*(1+prop);
    tau_w_ss = tau_w_ss*(1+prop);
    tau_K_ss = tau_K_ss*(1+prop);

    disp(['Solving for : ', S_.options.fiscal])
 
    %%%%%%%%%%%%%%%% Case gyss_b %%%%%%%%%%%%%%%%
elseif strcmp('gyss_tyss_b',S_.options.fiscal)

    [y_steady, Out_ ] = steadystate_rbc_gyss_tyss_b(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    by_ss     = log(y_steady(l_)); l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.ty_ss*y_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tyss_tauc(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tau_c_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.ty_ss*y_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tyss_tauw(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tau_w_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.ty_ss*y_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tyss_tauK(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    tau_K_ss  = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.ty_ss*y_ss;

    disp(['Solving for : ', S_.options.fiscal])

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

    [y_steady, Out_ ] = steadystate_rbc_gyss_tyss_taxprop(S_);

    if Out_.iterations > max_iter_steady || isempty(Out_.iterations)
        disp('Error in finding steady state')
        return
    end

    y_steady = exp(y_steady); l_ = 1;

    c_ss      = y_steady(l_);      l_ = l_ + 1;
    invest_ss = y_steady(l_);      l_ = l_ + 1;
    lambda_ss = y_steady(l_);      l_ = l_ + 1;
    r_F_ss    = y_steady(l_);      l_ = l_ + 1;
    kk_ss     = y_steady(l_);      l_ = l_ + 1;
    l_ss      = y_steady(l_);      l_ = l_ + 1;
    w_ss      = y_steady(l_);      l_ = l_ + 1;
    r_K_ss    = y_steady(l_);      l_ = l_ + 1;
    y_ss      = y_steady(l_);      l_ = l_ + 1;
    nx_ss     = log(y_steady(l_)); l_ = l_ + 1;
    ca_ss     = log(y_steady(l_)); l_ = l_ + 1;
    r_star_ss = y_steady(l_);      l_ = l_ + 1;
    b_F_ss    = log(y_steady(l_)); l_ = l_ + 1;
    r_ss      = y_steady(l_);      l_ = l_ + 1;
    b_ss      = log(y_steady(l_)); l_ = l_ + 1;
    g_ss      = y_steady(l_);      l_ = l_ + 1;
    prop      = y_steady(l_);      l_ = l_ + 1;
    p_b_ss    = y_steady(l_);
    phi_ss    = lambda_ss;
    tax_ss    = gy_ss*y_ss + ((1+kappa_b*p_b_ss)/z_ss - p_b_ss)*b_ss;

    tr_ss = S_.param.ty_ss*y_ss;

    tau_c_ss = tau_c_ss*(1-prop);
    tau_w_ss = tau_w_ss*(1-prop);
    tau_K_ss = tau_K_ss*(1-prop);

    disp(['Solving for : ', S_.options.fiscal])

end


%% Variable and dimensions

v_ = var_dim ; 

%% Structural matrices 

GAMMA0t  = zeros(v_.dim.n1+v_.dim.n2,v_.dim.n) ;
GAMMA1t  = zeros(v_.dim.n1+v_.dim.n2,v_.dim.n) ;
PSIt     = zeros(v_.dim.n1+v_.dim.n2,v_.dim.ns) ;
Ctilde   = zeros(v_.dim.n1+v_.dim.n2,1) ;

OMEGA    = diag([sig_z,sig_rstar,sig_zeta,sig_L,sig_I, sig_b, sig_g,sig_tauc,sig_tauw,sig_tauK,sig_tau]).^2;


% Equations

eqn = 1;

% Equation 1
% (z_ss-h*beta)(z_ss-h)lambda +  (z_ss-h*beta)(z_ss-h)/(1+tau_c_ss)*tau_c + (z_ss^2 + h^2*beta) c - h*beta*z_ss c(+1) + h*z_ss z
% - h*z*beta z(+1) - (z_ss-h)*z_ss zeta + (z_ss-h)*h*beta*rho_zeta zeta = h*z_ss c(-1) 

GAMMA0t(eqn,v_.lambda)= ((z_ss^sigma)-h*beta)*(z_ss-h);
GAMMA0t(eqn,v_.c)     = (sigma*z_ss^(sigma+1) + sigma*h^2*beta);
GAMMA0t(eqn,v_.E_c)   = - h*beta*sigma*z_ss ;
GAMMA1t(eqn,v_.c)     = h*sigma*z_ss^sigma;
GAMMA0t(eqn,v_.tau_c) = ((z_ss^sigma)-h*beta)*(z_ss-h)/(1+tau_c_ss);
GAMMA0t(eqn,v_.zeta)  = -(z_ss-h)*z_ss^sigma + (z_ss-h)*h*beta*rho_zeta;
GAMMA0t(eqn,v_.z)  = sigma*z_ss^(sigma)*h - sigma*h*z_ss*beta*rho_z;
Ctilde(eqn,1)         = ((z_ss^sigma)-h*beta)*(z_ss-h)*log(lambda_ss) + ((sigma*z_ss^(sigma+1) + sigma*h^2*beta)-(h*beta*sigma*z_ss)-h*sigma*z_ss^sigma)*log(c_ss) + ((z_ss^sigma)-h*beta)*(z_ss-h)/(1+tau_c_ss)*tau_c_ss + (sigma*z_ss^sigma*h - sigma*h*z_ss*beta*rho_z)*log(z_ss);

eqn = eqn + 1; 

% Equation 2
% lambda = lambda(+1) + r - rho_z z;

GAMMA0t(eqn,v_.lambda)    = 1;
GAMMA0t(eqn,v_.E_lambda)  = -1;
GAMMA0t(eqn,v_.r)         = -1;
GAMMA0t(eqn,v_.z)         = sigma*rho_z;
Ctilde(eqn,1)             = -log(r_ss) + sigma*rho_z*log(z_ss);

eqn = eqn + 1;

% Equation 2a
% 0 = p_b - kappa_b/(1+kappa_b*p_b_ss)*p_b(+1) + r
GAMMA0t(eqn,v_.p_b)      = 1;
GAMMA0t(eqn,v_.E_p_b)    = - kappa_b/(1+kappa_b*p_b_ss);
GAMMA0t(eqn,v_.r)        = 1;
Ctilde(eqn,1)            = log(p_b_ss) - kappa_b/(1+kappa_b*p_b_ss)*log(p_b_ss) + log(r_ss);

eqn = eqn + 1;


% Equation 3
% lambda = lambda(+1) + r_F - rho_z z;

GAMMA0t(eqn,v_.lambda)   = 1;
GAMMA0t(eqn,v_.E_lambda) = -1;
GAMMA0t(eqn,v_.r_F)      = -1;
GAMMA0t(eqn,v_.z)        = sigma*rho_z;
Ctilde(eqn,1)            = -log(r_F_ss) + sigma*rho_z*log(z_ss);

eqn = eqn + 1;

% Equation 4
% -phi + beta(1-tau_K_ss)/z_ss*r_K_ss*( lambda(+1) - rho_z z + r_K(+1)) - (beta r_K_ss/z_ss) tau_K(+1) + beta/z_ss*(1-delta)*(phi(+1) - rho_z z) = 0;
GAMMA0t(eqn,v_.phi)      = -1;
GAMMA0t(eqn,v_.E_lambda) = beta*(1-tau_K_ss)*r_K_ss/(z_ss^sigma);
GAMMA0t(eqn,v_.z)        = -sigma*rho_z*beta*(1-tau_K_ss)*r_K_ss/(z_ss^sigma) - sigma*rho_z*beta/(z_ss^sigma)*(1-delta) ;
GAMMA0t(eqn,v_.E_r_K)    = beta*(1-tau_K_ss)*r_K_ss/(z_ss^sigma);
GAMMA0t(eqn,v_.E_tau_K)  = -beta*r_K_ss/(z_ss^sigma);
GAMMA0t(eqn,v_.E_phi)    = beta*(1-delta)/(z_ss^sigma);
Ctilde(eqn,1)            = -log(phi_ss) + beta*(1-tau_K_ss)*r_K_ss/(z_ss^sigma)*log(lambda_ss) + (-sigma*rho_z*beta*(1-tau_K_ss)*r_K_ss/(z_ss^sigma) - sigma*rho_z*beta/(z_ss^sigma)*(1-delta))*log(z_ss)...
                         + (beta*(1-tau_K_ss)*r_K_ss/(z_ss^sigma))*log(r_K_ss) -(beta*r_K_ss/(z_ss^sigma))*tau_K_ss + beta/(z_ss^sigma)*(1-delta)*log(phi_ss);

eqn = eqn + 1;


% Equation 5
% -lambda + phi - z_ss^2 * upsilon*z - z_ss^2*upsilon*(1 + bet)*invest = -z_ss^2*upsilon*(invest(-1) + beta*invest(+1) + beta*z(+1))

GAMMA0t(eqn,v_.lambda)   = -1;
GAMMA0t(eqn,v_.phi)      = 1;
GAMMA0t(eqn,v_.invest)   = -z_ss^2*Upsilon*(1 + beta*z_ss^(1-sigma));
GAMMA0t(eqn,v_.z)        = (beta*rho_z*z_ss^(1-sigma)-1)*z_ss^2*Upsilon;
GAMMA1t(eqn,v_.invest)   = -z_ss^2*Upsilon;
GAMMA0t(eqn,v_.E_invest) = z_ss^2*Upsilon*beta*z_ss^(1-sigma);
GAMMA0t(eqn,v_.zeta_I)   = 1;
Ctilde(eqn,1)            = (beta*rho_z*z_ss^(1-sigma)-1)*z_ss^2*Upsilon*log(z_ss);

eqn = eqn + 1;

% Equation 6
% -lambda = -1/(1-tau_w_ss)*tau_w + w - nu*l - zeta(t) - zeta_L(t)

GAMMA0t(eqn,v_.lambda) = -1;
GAMMA0t(eqn,v_.tau_w)  = 1/(1-tau_w_ss);
GAMMA0t(eqn,v_.w)      = -1;
GAMMA0t(eqn,v_.l)      = nu;
GAMMA0t(eqn,v_.zeta)   = 1;
GAMMA0t(eqn,v_.zeta_L) = 1;
Ctilde(eqn,1)       = -log(lambda_ss) + 1/(1-tau_w_ss)*(tau_w_ss) - log(w_ss) + nu*log(l_ss) ;

eqn = eqn + 1;

% Equation 7
% kk  = (1-delta)(z_ss^-1)*kk(-1) + invest_ss/kk_ss*invest - (1-delta)(z_ss^-1)*z;

GAMMA0t(eqn,v_.kk)     = 1;
GAMMA1t(eqn,v_.kk)     = (1-delta)*(z_ss^-1);
GAMMA0t(eqn,v_.z)      = (1-delta)*(z_ss^-1) ;
GAMMA0t(eqn,v_.invest) = - (z_ss - 1 + delta)/z_ss;
GAMMA0t(eqn,v_.zeta_I) = - (z_ss - 1 + delta)/z_ss; 
Ctilde(eqn,1)          = log(kk_ss) - (1-delta)*(z_ss^-1)*log(kk_ss) - ((z_ss - 1 + delta)/z_ss)*log(invest_ss) + (1-delta)*(z_ss^-1)*log(z_ss);

eqn = eqn + 1; 

% Equation 8
% y = a + alpha*k(-1) + (1-alpha)*l - alpha*z;

GAMMA0t(eqn,v_.y)  = 1;
% GAMMA0t(eqn,v_.a)  = -1;
GAMMA1t(eqn,v_.kk) = alpha;
GAMMA0t(eqn,v_.l)  = -(1-alpha);
GAMMA0t(eqn,v_.z)  = alpha;
Ctilde(eqn,1)      = log(y_ss) - log(A_ss) - alpha*log(kk_ss) - (1-alpha)*log(l_ss) + alpha*log(z_ss);

eqn = eqn + 1;

% Equation 9
% r_K - y - z = - kk(-1);
GAMMA0t(eqn,v_.r_K)  = 1;
GAMMA0t(eqn,v_.y)    = -1;
GAMMA1t(eqn,v_.kk)   = -1;
GAMMA0t(eqn,v_.z)    = -1;
Ctilde(eqn,1)        = log(r_K_ss) - log(y_ss) + log(kk_ss) - log(z_ss);

eqn = eqn + 1;

% Equation 10
% y = w + l;
GAMMA0t(eqn,v_.y)   = 1;
GAMMA0t(eqn,v_.w)   = -1;
GAMMA0t(eqn,v_.l)   = -1;
Ctilde(eqn,1)       = log(y_ss) - log(w_ss) - log(l_ss) ;

eqn = eqn + 1;

% Equation 11
% nx = y_ss*y - c_ss*c - i_ss*i - g_ss*g ;
GAMMA0t(eqn,v_.nx)     = 1;
GAMMA0t(eqn,v_.y)      = -y_ss;
GAMMA0t(eqn,v_.c)      = c_ss;
GAMMA0t(eqn,v_.invest) = invest_ss;
GAMMA0t(eqn,v_.g)      = g_ss;
Ctilde(eqn,1)          = nx_ss - y_ss*log(y_ss) + c_ss*log(c_ss) + invest_ss*log(invest_ss) + g_ss*log(g_ss);

eqn = eqn + 1;

% Equation 12
% ca = nx + r_F_ss*b_F_ss/z_ss r_F(-1) + (r_F_ss-1)/z_ss b_F(-1) - (r_F_ss-1)/z_ss z;
GAMMA0t(eqn,v_.ca)   = 1;
GAMMA0t(eqn,v_.nx)   = -1;
GAMMA1t(eqn,v_.b_F)  = (r_F_ss-1)/z_ss;
GAMMA1t(eqn,v_.r_F)  = r_F_ss*b_F_ss/z_ss;
Ctilde(eqn,1)        = ca_ss - nx_ss - r_F_ss*b_F_ss/z_ss * log(r_F_ss) - (r_F_ss-1)/z_ss*b_F_ss  + (r_F_ss-1)/z_ss *log(z_ss);

eqn = eqn + 1;

% Equation 13
% b_F = nx + (r_F_ss/z_ss) b_F(-1) + (r_F_ss b_F_ss/z_ss) r_F(-1) - (r_F_ss b_F_ss/z_ss) z;
GAMMA0t(eqn,v_.b_F)  = 1;
GAMMA0t(eqn,v_.nx)   = -1;
GAMMA0t(eqn,v_.z)    = r_F_ss*b_F_ss/z_ss ;
GAMMA1t(eqn,v_.b_F)  = r_F_ss/z_ss;
GAMMA1t(eqn,v_.r_F)  = r_F_ss*b_F_ss/z_ss;
Ctilde(eqn,1)        = (1-r_F_ss/z_ss)*b_F_ss - nx_ss - r_F_ss*b_F_ss/z_ss * (log(r_F_ss)-log(z_ss)); 

eqn = eqn + 1;

% Equation 14
% r_F - r_star -zeta_b + psi_b/y_ss*b_F - psi_b b_F_ss/y_ss y = 0;
GAMMA0t(eqn,v_.r_F)    = 1;
GAMMA0t(eqn,v_.r_star) = -1;
GAMMA0t(eqn,v_.zeta_b) = -1;
GAMMA0t(eqn,v_.b_F)    = psi_b/y_ss;
GAMMA0t(eqn,v_.y)      = - psi_b*b_F_ss/y_ss ; 
Ctilde(eqn,1)          = log(r_F_ss) - log(r_star_ss) - psi_b*b_F_ss/y_ss * log(y_ss) + psi_b/y_ss*b_F_ss - zeta_b_ss;

eqn = eqn + 1;

% Equation 15
% g_ss * g  -(1+kappa_b*p_b_ss)*b_ss/z_ss * z + (1+kappa_b*p_b_ss)/z_ss b(-1) + (kappa_b*p_b_ss/z_ss*b_ss - p_b_ss*b_ss)*p_b =  b + c_ss tau_c
% tau_c_ss*c_ss * c + w_ss l_ss tau_w + tau_w_ss*w_ss*l_ss* (w + l) 
% r_K_ss k_ss/z_ss tau_K + (tau_K_ss r_K_ss k_ss/z_ss) (r_K + k(-1) - z) + tau;

GAMMA0t(eqn,v_.g)      = g_ss;
GAMMA1t(eqn,v_.b)      = -(1+kappa_b*p_b_ss)/z_ss;
GAMMA0t(eqn,v_.z)      = -(1+kappa_b*p_b_ss)*b_ss/z_ss + (tau_K_ss*r_K_ss*kk_ss/z_ss);
GAMMA0t(eqn,v_.b)      = -p_b_ss;
GAMMA0t(eqn,v_.p_b)    = kappa_b*p_b_ss*b_ss/z_ss - p_b_ss*b_ss;
GAMMA0t(eqn,v_.tau_c)  = -c_ss;
GAMMA0t(eqn,v_.c)      = -tau_c_ss*c_ss;
GAMMA0t(eqn,v_.tau_w)  = -w_ss*l_ss;
GAMMA0t(eqn,v_.w)      = -tau_w_ss*w_ss*l_ss;
GAMMA0t(eqn,v_.l)      = -tau_w_ss*w_ss*l_ss;
GAMMA0t(eqn,v_.tau_K)  = -r_K_ss*kk_ss/z_ss;
GAMMA0t(eqn,v_.r_K)    = -tau_K_ss*r_K_ss*kk_ss/z_ss;
GAMMA1t(eqn,v_.kk)     = tau_K_ss*r_K_ss*kk_ss/z_ss;
GAMMA0t(eqn,v_.tr)     = -1;
Ctilde(eqn,1)          = g_ss*log(g_ss) + (-(1+kappa_b*p_b_ss)*b_ss/z_ss + (tau_K_ss*r_K_ss*kk_ss/z_ss))*log(z_ss) + ((1+kappa_b*p_b_ss)/z_ss -p_b_ss) * b_ss...
                        -c_ss*tau_c_ss - tau_c_ss*c_ss * log(c_ss) - w_ss*l_ss*tau_w_ss - tau_w_ss*w_ss*l_ss * (log(w_ss) + log(l_ss)) ... 
                        - (r_K_ss*kk_ss/z_ss)*tau_K_ss - tau_K_ss*r_K_ss*kk_ss/z_ss * (log(r_K_ss)+log(kk_ss)) - tr_ss + (kappa_b*p_b_ss*b_ss/z_ss - p_b_ss*b_ss)*log(p_b_ss) ;

eqn = eqn + 1;

% Equation 16 (defines total tax revenue)
% tax_ss * tax = c_ss tau_c
% tau_c_ss*c_ss * c + w_ss l_ss tau_w + tau_w_ss*w_ss*l_ss* (w + l) 
% r_K_ss k_ss/z_ss tau_K + (tau_K_ss r_K_ss k_ss/z_ss) (r_K + k(-1) - z) + tau;
GAMMA0t(eqn,v_.tax)    =  tax_ss;
GAMMA0t(eqn,v_.z)      = (tau_K_ss*r_K_ss*kk_ss/z_ss);
GAMMA0t(eqn,v_.tau_c)  = -c_ss;
GAMMA0t(eqn,v_.c)      = -tau_c_ss*c_ss;
GAMMA0t(eqn,v_.tau_w)  = -w_ss*l_ss;
GAMMA0t(eqn,v_.w)      = -tau_w_ss*w_ss*l_ss;
GAMMA0t(eqn,v_.l)      = -tau_w_ss*w_ss*l_ss;
GAMMA0t(eqn,v_.tau_K)  = -r_K_ss*kk_ss/z_ss;
GAMMA0t(eqn,v_.r_K)    = -tau_K_ss*r_K_ss*kk_ss/z_ss;
GAMMA1t(eqn,v_.kk)     = tau_K_ss*r_K_ss*kk_ss/z_ss;
GAMMA0t(eqn,v_.tr)     = -1;
Ctilde(eqn,1)          = tax_ss*log(tax_ss) + ((tau_K_ss*r_K_ss*kk_ss/z_ss))*log(z_ss) ...
                         -c_ss*tau_c_ss - tau_c_ss*c_ss * log(c_ss) - w_ss*l_ss*tau_w_ss - tau_w_ss*w_ss*l_ss * (log(w_ss) + log(l_ss)) ... 
                         - (r_K_ss*kk_ss/z_ss)*tau_K_ss - tau_K_ss*r_K_ss*kk_ss/z_ss * (log(r_K_ss)+log(kk_ss)) - tr_ss;

eqn = eqn + 1;

% Equation 16
% g = (1-rho_g)*log(g_ss) - (1-rho_g)*gam_gb*b_ss/y_ss*log(y_ss) + rho_g*g(-1) - (1 - rho_g)*(gam_gb/y_ss)*b(-1) + (1 - rho_g)*(gam_gb*b_ss/y_ss) * y(-1)) + (1-rho_g)*gam_gb by_star +
%     - (1-rho_g)*gam_gy * ((y(-1) - log(y_ss)) +  vg;

GAMMA0t(eqn,v_.g)       = 1;
GAMMA1t(eqn,v_.y)       = (1 - rho_g- rho2_g)*(gam_gb*p_b_ss*b_ss/y_ss - gam_gy);
GAMMA1t(eqn,v_.p_b)     = - (1 - rho_g- rho2_g)*(gam_gb*p_b_ss*b_ss/y_ss - gam_gy);
GAMMA1t(eqn,v_.g)       = rho_g;
GAMMA1t(eqn,v_.lag1_g)  = rho2_g;
GAMMA1t(eqn,v_.b)       = -(1 - rho_g - rho2_g)*gam_gb*p_b_ss/y_ss;
GAMMA0t(eqn,v_.vg)      = -1;
GAMMA0t(eqn,v_.by_star) = -(1 - rho_g - rho2_g)*gam_gb ; 
Ctilde(eqn)             = (1 - rho_g - rho2_g)*log(g_ss) - (1 - rho_g - rho2_g)*(gam_gb*p_b_ss*b_ss/y_ss-gam_gy)*(log(y_ss) - log(p_b_ss));

eqn = eqn + 1;

% 
GAMMA0t(eqn,v_.lag1_g)      = 1;
GAMMA1t(eqn,v_.g)           = 1;

eqn = eqn + 1;

% Equation 17
% tau_c = (1-rho_c)*tau_c_ss + (1-rho_c)*gam_cb*b_ss/y_ss*log(y_ss) + rho_c*tau_c(-1)  + (1 - rho_c)*gam_cb/y_ss*b(-1) - (1 - rho_c)*gam_cb*b_ss/y_ss * y(-1)) - (1-rho_c)*gam_cb by_star +  ec;
%       + (1-rho_c)*gam_cy * ((y(-1) - log(y_ss))

GAMMA0t(eqn,v_.tau_c)     = 1;
GAMMA1t(eqn,v_.y)         = -(1 - rho_c- rho2_c)*(gam_cb*p_b_ss*b_ss/y_ss - gam_cy);
GAMMA1t(eqn,v_.p_b)        = (1 - rho_c- rho2_c)*(gam_cb*b_ss*p_b_ss/y_ss - gam_cy);
GAMMA1t(eqn,v_.tau_c)     = rho_c;
GAMMA1t(eqn,v_.lag1_tau_c)= rho2_c;
GAMMA1t(eqn,v_.b)         = (1 - rho_c- rho2_c)*gam_cb*p_b_ss/y_ss;
GAMMA0t(eqn,v_.by_star)   = (1 - rho_c- rho2_c)*gam_cb ; 
PSIt(eqn,v_.shock.ectau)  = 1;
Ctilde(eqn)               = (1 - rho_c- rho2_c)*tau_c_ss + (1 - rho_c- rho2_c)*(gam_cb*p_b_ss*b_ss/y_ss-gam_cy)*(log(y_ss) - log(p_b_ss));

eqn = eqn + 1;

% 
GAMMA0t(eqn,v_.lag1_tau_c)      = 1;
GAMMA1t(eqn,v_.tau_c)           = 1;

eqn = eqn + 1;

% Equation 18
% tau_w  = (1-rho_w)*(tau_w_ss) + (1-rho_w)*gam_wb*b_ss/y_ss*log(y_ss) + rho_w*tau_w(-1) + (1 - rho_w)*gam_wb/y_ss*b(-1) - (1 - rho_w)*gam_wb*b_ss/y_ss * y(-1)) - (1-rho_w)*gam_wb by_star + ew;
%       + (1-rho_w)*gam_wy * ((y(-1) - log(y_ss))

GAMMA0t(eqn,v_.tau_w)    = 1;
GAMMA1t(eqn,v_.y)        = -(1 - rho_w - rho2_w)*(gam_wb*b_ss*p_b_ss/y_ss - gam_wy);
GAMMA1t(eqn,v_.p_b)     = (1 - rho_w - rho2_w)*(gam_wb*p_b_ss*b_ss/y_ss - gam_wy);
GAMMA1t(eqn,v_.tau_w)    = rho_w;
GAMMA1t(eqn,v_.lag1_tau_w)= rho2_w;
GAMMA1t(eqn,v_.b)        = (1 - rho_w - rho2_w)*gam_wb*p_b_ss/y_ss;
GAMMA0t(eqn,v_.by_star)  = (1 - rho_w - rho2_w)*gam_wb ; 
PSIt(eqn,v_.shock.ewtau) = 1;
Ctilde(eqn)              = (1 - rho_w - rho2_w)*(tau_w_ss) + (1 - rho_w - rho2_w)*(gam_wb*p_b_ss*b_ss/y_ss - gam_wy)*(log(y_ss) - log(p_b_ss));

eqn = eqn + 1;

% 
GAMMA0t(eqn,v_.lag1_tau_w)      = 1;
GAMMA1t(eqn,v_.tau_w)           = 1;

eqn = eqn + 1;

% Equation 19
% tau_K  = constant + rho_K*tau_K(-1) + (1 - rho_K)*gam_Kb/y_ss*b(-1) - (1 - rho_K)*gam_Kb*b_ss/y_ss * y(-1)) - (1-rho_K)*gam_Kb by_star + eK;
%       + (1-rho_K)*gam_Ky * ((y(-1) - log(y_ss))

GAMMA0t(eqn,v_.tau_K)    = 1;
GAMMA1t(eqn,v_.y)        = -(1 - rho_K - rho2_K)*(gam_Kb*b_ss*p_b_ss/y_ss - gam_Ky);
GAMMA1t(eqn,v_.p_b)     = (1 - rho_K - rho2_K)*(gam_Kb*p_b_ss*b_ss/y_ss - gam_Ky);
GAMMA1t(eqn,v_.tau_K)    = rho_K;
GAMMA1t(eqn,v_.lag1_tau_K)= rho2_K;
GAMMA1t(eqn,v_.b)        = (1 - rho_K - rho2_K)*gam_Kb*p_b_ss/y_ss;
GAMMA0t(eqn,v_.by_star)  = (1 - rho_K - rho2_K)*gam_Kb ; 
PSIt(eqn,v_.shock.eKtau) = 1;
Ctilde(eqn)              = (1 - rho_K - rho2_K)*(tau_K_ss) + (1 - rho_K - rho2_K)*(gam_Kb*p_b_ss*b_ss/y_ss - gam_Ky)*(log(y_ss) - log(p_b_ss));

eqn = eqn + 1;

% 
GAMMA0t(eqn,v_.lag1_tau_K)      = 1;
GAMMA1t(eqn,v_.tau_K)           = 1;

eqn = eqn + 1;

% Equation 20
% tau  = const + rho_tau*tau(-1) + (1 - rho_tau)*gam_taub/y_ss*b(-1) - (1 - rho_tau)*gam_taub*b_ss/y_ss * y(-1)) - (1-rho_tau)*gam_taub by_star + etau;
%       + (1-rho_K)*gam_tauy * ((y(-1) - log(y_ss))

GAMMA0t(eqn,v_.tr)      = 1;
GAMMA1t(eqn,v_.y)       = -(1 - rho_tau -rho2_tau)*(gam_taub*p_b_ss*b_ss/y_ss - gam_tauy);
GAMMA1t(eqn,v_.p_b)     = (1 - rho_tau -rho2_tau)*(gam_taub*p_b_ss*b_ss/y_ss - gam_tauy);
GAMMA1t(eqn,v_.tr)      = rho_tau;
GAMMA1t(eqn,v_.lag1_tr)= rho2_tau;
GAMMA1t(eqn,v_.b)       = (1 - rho_tau -rho2_tau)*gam_taub*p_b_ss/y_ss;
GAMMA0t(eqn,v_.by_star) = (1 - rho_tau -rho2_tau)*gam_taub ; 
PSIt(eqn,v_.shock.etau) = 1;
Ctilde(eqn)             = (1 - rho_tau -rho2_tau)*(tr_ss) + (1 - rho_tau -rho2_tau)*(gam_taub*p_b_ss*b_ss/y_ss-gam_tauy)*(log(y_ss) - log(p_b_ss));

eqn = eqn + 1;

% 
GAMMA0t(eqn,v_.lag1_tr)      = 1;
GAMMA1t(eqn,v_.tr)           = 1;

eqn = eqn + 1;

% % Equation 21
% % a = rho*a(-1)+ea;
% 
% GAMMA0t(eqn,a)   = 1;
% GAMMA1t(eqn,a)   = rho_a;
% PSIt(eqn,ea)     = 1;
% Ctilde(eqn,1)    = (1-rho_a)*log(A_ss);
% 
% eqn = eqn + 1;

% Equation 22
% z = rho_z z(-1) + ez;

GAMMA0t(eqn,v_.z)     = 1;
GAMMA1t(eqn,v_.z)     = rho_z;
PSIt(eqn,v_.shock.ez) = 1;
Ctilde(eqn,1)         = (1-rho_z)*log(z_ss);

eqn = eqn + 1;

% Equation 23
% r_star = rho*r_star(-1) + erstar;

GAMMA0t(eqn,v_.r_star)    = 1;
GAMMA1t(eqn,v_.r_star)    = rho_r_star;
PSIt(eqn,v_.shock.erstar) = 1;
Ctilde(eqn,1)             = (1-rho_r_star)*log(r_star_ss);

eqn = eqn + 1;

% Equation 24
% zeta = rho_zeta zeta(-1) + ezeta;

GAMMA0t(eqn,v_.zeta)     = 1;
GAMMA1t(eqn,v_.zeta)     = rho_zeta;
PSIt(eqn,v_.shock.ezeta) = 1;

eqn = eqn + 1;

% Equation 25
% zeta_L = rho_L zeta_L(-1) + eL;

GAMMA0t(eqn,v_.zeta_L)  = 1;
GAMMA1t(eqn,v_.zeta_L)  = rho_L;
PSIt(eqn,v_.shock.eL)   = 1;

eqn = eqn + 1;

% Equation 26
% zeta_I = rho_I zeta_I(-1) + eI;

GAMMA0t(eqn,v_.zeta_I)  = 1;
GAMMA1t(eqn,v_.zeta_I)  = rho_I;
PSIt(eqn,v_.shock.eI)   = 1;

eqn = eqn + 1;

% Equation 27
% zeta_b = (1-rho_b) zeta_b_ss + rho_b zeta_b(-1) + eb;

GAMMA0t(eqn,v_.zeta_b) = 1;
GAMMA1t(eqn,v_.zeta_b) = rho_b;
PSIt(eqn,v_.shock.eb)      = 1;
Ctilde(eqn,1)          = (1-rho_b)*zeta_b_ss; 

eqn = eqn + 1;

% Equation 28
% vg = rho_vg vg(-1) + eg;

GAMMA0t(eqn,v_.vg)    = 1;
GAMMA1t(eqn,v_.vg)    = rho_vg;
PSIt(eqn,v_.shock.eg) = 1;

eqn = eqn + 1;

% Equation 29
% by_star = (1-rho_by) b_ss/y_ss + rho_by by_star(-1) + eby;

GAMMA0t(eqn,v_.by_star) = 1;
GAMMA1t(eqn,v_.by_star) = rho_by;
% PSIt(eqn,v_.shock.eby)  = 1;
Ctilde(eqn,1)           = (1-rho_by)*p_b_ss*b_ss/y_ss;

eqn = eqn + 1;

% Observation equations

% Equation 30 dy_obs
GAMMA0t(eqn,v_.dy_obs) =  1;
GAMMA0t(eqn,v_.y)      = -1;
GAMMA1t(eqn,v_.y)      = -1;
GAMMA0t(eqn,v_.z)      = -1;

eqn = eqn + 1; 

% Equation 31 dinvest_obs
GAMMA0t(eqn,v_.dinvest_obs) = 1;
GAMMA0t(eqn,v_.invest)      = -1;
GAMMA1t(eqn,v_.invest)      = -1;
GAMMA0t(eqn,v_.z)           = -1;

eqn = eqn + 1; 

% Equation 31 dinvest_obs
GAMMA0t(eqn,v_.dc_obs) = 1;
GAMMA0t(eqn,v_.c)      = -1;
GAMMA1t(eqn,v_.c)      = -1;
GAMMA0t(eqn,v_.z)      = -1;

eqn = eqn + 1; 

% Equation 32 nx_obs is (nx to y ratio)
GAMMA0t(eqn,v_.nx_obs) = 1; 
GAMMA0t(eqn,v_.nx)     = -1/y_ss;
GAMMA0t(eqn,v_.y)      = nx_ss/y_ss;
Ctilde(eqn,1)          = (nx_ss/y_ss) * log(y_ss);   

eqn = eqn + 1; 

% Equation 33 dw_obs
GAMMA0t(eqn,v_.dw_obs) = 1; 
GAMMA0t(eqn,v_.w)      = -1;
GAMMA1t(eqn,v_.w)      = -1;
GAMMA0t(eqn,v_.z)      = -1;

eqn = eqn + 1; 

% Equation 34 r_obs
GAMMA0t(eqn,v_.r_obs) = 1; 
GAMMA0t(eqn,v_.r)     = -1;

eqn = eqn + 1; 

% Equation 35 r_star_obs
GAMMA0t(eqn,v_.r_star_obs) = 1; 
GAMMA0t(eqn,v_.r_star)     = -1;

eqn = eqn + 1; 

% Equation 36 gy_obs
GAMMA0t(eqn,v_.gy_obs) = 1;
GAMMA0t(eqn,v_.g)      = -1 ; 
GAMMA0t(eqn,v_.y)      = 1 ; 
Ctilde(eqn,1)          = 0;

eqn = eqn + 1; 

% Equation 37 by_obs
GAMMA0t(eqn,v_.by_obs) = 1;
GAMMA0t(eqn,v_.p_b)    = - by_ss ;
GAMMA0t(eqn,v_.b)      = -p_b_ss/y_ss;
GAMMA0t(eqn,v_.y)      = by_ss ;
Ctilde(eqn,1)          = by_ss*log(y_ss) - by_ss*log(p_b_ss);

eqn = eqn + 1; 

% Equation 38 rev_c_y_obs
GAMMA0t(eqn,v_.rev_c_y_obs) = 1;
GAMMA0t(eqn,v_.tau_c)       = -c_ss/y_ss;
GAMMA0t(eqn,v_.c)           = -tau_c_ss*c_ss/y_ss;
GAMMA0t(eqn,v_.y)           = tau_c_ss*c_ss/y_ss;
Ctilde(eqn,1)               = (tau_c_ss*c_ss/y_ss) * (log(y_ss) - log(c_ss));

eqn = eqn + 1; 

% Equation 39 rev_w_y_obs
GAMMA0t(eqn,v_.rev_w_y_obs) = 1;
GAMMA0t(eqn,v_.tau_w)       = -w_ss*l_ss/y_ss;
GAMMA0t(eqn,v_.w)           = -tau_w_ss*w_ss*l_ss/y_ss;
GAMMA0t(eqn,v_.l)           = -tau_w_ss*w_ss*l_ss/y_ss;
GAMMA0t(eqn,v_.y)           = tau_w_ss*w_ss*l_ss/y_ss;
Ctilde(eqn,1)               = (tau_w_ss*w_ss*l_ss/y_ss)*(log(y_ss) - log(l_ss) - log(w_ss));

eqn = eqn + 1; 

% Equation 40 rev_K_y_obs
GAMMA0t(eqn,v_.rev_K_y_obs) = 1;
GAMMA0t(eqn,v_.tau_K)       = -r_K_ss*kk_ss/(z_ss*y_ss);
GAMMA0t(eqn,v_.r_K)         = -tau_K_ss*r_K_ss*kk_ss/(z_ss*y_ss);
GAMMA1t(eqn,v_.kk)          = tau_K_ss*r_K_ss*kk_ss/(z_ss*y_ss);
GAMMA0t(eqn,v_.y)           = tau_K_ss*r_K_ss*kk_ss/(z_ss*y_ss);
GAMMA0t(eqn,v_.z)           = tau_K_ss*r_K_ss*kk_ss/(z_ss*y_ss);
Ctilde(eqn,1)               = (tau_K_ss*r_K_ss*kk_ss/(z_ss*y_ss))*(log(y_ss)+log(z_ss)-log(kk_ss)-log(r_K_ss));

eqn = eqn + 1; 

% Equation 41 s
GAMMA0t(eqn,v_.s) = 1;
GAMMA0t(eqn,v_.g) = gy_ss*y_ss;
GAMMA0t(eqn,v_.tax) = -tax_ss;
Ctilde(eqn,1)  = (tax_ss - g_ss) + gy_ss*y_ss*log(g_ss) - tax_ss*log(tax_ss);

%% Structural matrics for Smats.m
% Add in Expectational Errors

GAM0    = [ GAMMA0t; [zeros(v_.dim.k,v_.dim.n1), eye(v_.dim.k), zeros(v_.dim.k,v_.dim.n2)] ];
GAM1    = [ GAMMA1t; [zeros(v_.dim.k,v_.dim.n1+v_.dim.n2),eye(v_.dim.k)] ];
PSI     = [ PSIt; zeros(v_.dim.k,v_.dim.ns) ];
PPI     = [ zeros(v_.dim.n1+v_.dim.n2,v_.dim.k); eye(v_.dim.k) ];
C_       = [ Ctilde; zeros(v_.dim.k,1) ];

st_mat.GAM0 = GAM0;
st_mat.GAM1 = GAM1;
st_mat.PSI  = PSI;
st_mat.PPI  = PPI;
st_mat.C_   = C_;
st_mat.OMEGA = OMEGA; 
st_mat.y_steady = y_steady;

