% This function computes additional variables for the Home bias model

function td = compute_addl_vars_StoneGeary(P,td)

% tau
td.tau      = td.r-td.rstar;
a_ss     = P.h_inv(P.rho*P.b0);

td.x        = exp(td.a);
x_ss        = exp(a_ss);

td.x_ss = (td.x-x_ss)/x_ss;

td.p        = P.p_fn(td.x);
p_ss        = P.p_fn(x_ss);
td.p_ss     = (td.p - p_ss)/p_ss;

td.cT_h2m   = P.c_lbar + P.alpha * (td.p*P.yNT_h2m + P.yT_h2m - P.c_lbar);
cT_h2m_ss   = P.c_lbar + P.alpha * (p_ss*P.yNT_h2m + P.yT_h2m - P.c_lbar);
td.cNT_h2m  = (1-P.alpha) * (td.p*P.yNT_h2m + P.yT_h2m - P.c_lbar) ./ td.p;
cNT_h2m_ss  = (1-P.alpha) * (p_ss*P.yNT_h2m + P.yT_h2m - P.c_lbar) / p_ss;

td.cNT_R    = (1-P.alpha)/P.alpha * td.x / td.p;
cNT_R_ss    = (1-P.alpha)/P.alpha * x_ss / p_ss;
td.cT_R     = td.x + P.c_lbar;
cT_R_ss     = x_ss + P.c_lbar;

td.z_h2m = td.cT_h2m - P.c_lbar + td.p .* td.cNT_h2m;  % total expenditure h2m
z_h2m_ss = cT_h2m_ss - P.c_lbar + p_ss .* cNT_h2m_ss;
td.z_R   = td.cT_R - P.c_lbar + td.p .* td.cNT_R;  % total expenditure R
z_R_ss   = cT_R_ss - P.c_lbar + p_ss .* cNT_R_ss;

td.z_h2m_ss = (td.z_h2m - z_h2m_ss)/z_h2m_ss;
td.z_R_ss   = (td.z_R - z_R_ss)/z_R_ss;

td.CPI      = (1/P.alpha)^P.alpha * (td.p/(1-P.alpha)).^(1-P.alpha);
td.RER      = P.RER_fn(td.a);

RER_ss   = P.RER_fn(a_ss); % original ss
td.RERss = td.RER / RER_ss - 1;

td.y_dollar = P.yT + td.p * P.yNT;
td.y_PPP    = td.y_dollar./td.CPI;

td.reserves = td.b + td.bI + td.bIH;
td.nfa      = td.b;

% break asymmetry if tau = 0 ?
if strcmp(P.break_indeterminacy,'res_const') || strcmp(P.break_indeterminacy,'res_linear')
    
    % go over stretches of tau = 0
    ind_tau_0 = find( abs(td.tau) < 1e-10);
    nonzero_tau_logical = abs(td.tau) >= 1e-10;
    
    for i = ind_tau_0
        last_i_nonzero_tau = find( nonzero_tau_logical & ([1:length(td.t)]<i),1,'last');
        next_i_nonzero_tau = find( nonzero_tau_logical & ([1:length(td.t)]>i),1,'first');
        if isempty(last_i_nonzero_tau)
            last_i_nonzero_tau = 1;
        end
        if isempty(next_i_nonzero_tau)
            next_i_nonzero_tau = length(td.t);
        end
        
        if strcmp(P.break_indeterminacy,'res_const')
            td.reserves_new(i) = td.reserves(last_i_nonzero_tau);
        end
        
        if strcmp(P.break_indeterminacy,'res_linear')
            weight = (td.t(i) - td.t(last_i_nonzero_tau)) / (td.t(next_i_nonzero_tau) - td.t(last_i_nonzero_tau));
            td.reserves_new(i) = (1-weight)*td.reserves(last_i_nonzero_tau) + weight*td.reserves(next_i_nonzero_tau);
        end
    end
    
    % adjust bI and bIH
    % for simplicity: adjust proportionately as to how much "slack" i.e. unused masses of GAmma = 0 they
    % still have
    for i = ind_tau_0
        res_diff = td.reserves_new(i)- td.reserves(i);
        if res_diff ~= 0
            if res_diff > 0
                % intermediaries hold more stuff
                small_r = 1e-9;
            elseif res_diff < 0
                % intermediaries hold less stuff
                small_r = -1e-9;
            end
            slack_F = P.D.fn(td.rstar(i)+small_r,td.rstar(i),td.q(i),td.qstar(i),td.t(i));
            slack_H = P.D.Hfn(td.rstar(i)+small_r,td.rstar(i),td.q(i),td.qstar(i),td.t(i));
            weight_F = slack_F / (slack_H + slack_F);

            td.bI(i)  = td.bI(i) + weight_F * res_diff;
            td.bIH(i) = td.bIH(i) + (1-weight_F) * res_diff;
            td.reserves(i) = td.reserves_new(i);
            assert( abs(- td.reserves(i) + td.b(i) + td.bI(i) + td.bIH(i))<1e-10)
        end
    end
end

% update objects
td.resGDP   = td.reserves/td.y_PPP(end); % reserves over initial GDP
td.privnfa = td.nfa - td.reserves;
td.privnfaGDP = td.privnfa/td.y_PPP(end); 
td.nfaGDP     = td.nfa/td.y_PPP(end);

% compute cost
td.cost     = td.tau.*td.bI;
if isfield(P.D,'xi')
    td.xi       = P.D.xi(td.t);
    td.tau_xi   = td.tau + td.xi/2;
end
td.costGDP  = td.cost/td.y_PPP(end); % relative to initial GDP
end