

% Find steady state in the safe model (sigma = sigma_lower_bar)
model_safe_paramfile_spreads
%param_core
%param_core_alt_large_tau_3;

%gammap=0.01005;
gammap = 0.03;

% define steady-state values for:
% c, lambda_c, d, r_d, r_e, r_l, w, y, h, k, invest, l, e, a, hr, kr, yr, r_lr, Tax

%\sigma_{0}d^{-\sigma_{d}}=c^{-\sigma_{c}}\left(1-\beta R^{d}\right)

% 
r_e_upper_bar = (1+kappab)/betap;

excess_risk_switch = 1;

sigma = sigma_bar;

algo = char('active-set', 'trust-region-reflective', 'interior-point', 'interior-point-convex', 'levenberg-marquardt', 'trust-region-dogleg', 'lm-line-search','sqp');
options = optimset('Display','none','Jacobian','off','MaxFunEvals',1e4,'MaxIter',1e4,'TolFun',1e-9,'TolX',1e-9,'Algorithm',algo(5,:));


%resid = function_r_l(r_e,betap,xip,deltap,alphap,sigma_lower_bar,gammap,taup,sigmap_zero,sigmap_d, sigmap_c,kappap,pip)

%changes_names = char('gammap');
%changes_vec = [gammap];
Tax=0;
Tax1=1;
%guess=[r_e;gammap;sigmap_zero];

guess=[r_e_upper_bar+0.001];

while abs(Tax-Tax1)>1e-9
    
    Tax1=Tax;
    
%     x = fsolve(@(x) function_r_l_fixed_spread(x,rl_rd_spreadp,betap,xip,deltap,alphap,sigma,sigma_lower_bar,sigma_bar,taup,sigmap_zero,sigmap_d, sigmap_c,kappap,pip,f_costs,kappab,Tax1),guess,options);
%     
%     [resid1] = function_r_l_fixed_spread(x,rl_rd_spreadp,betap,xip,deltap,alphap,sigma,sigma_lower_bar,sigma_bar,taup,sigmap_zero,sigmap_d, sigmap_c,kappap,pip,f_costs,kappab,Tax1);
%     
%     r_l=x(1);
%     gammap=x(2);
%     sigmap_zero = x(3);
    
    
    x = fsolve(@(x) function_r_l_risky(x,betap,xip,deltap,alphap,sigma,sigma_lower_bar,taup,sigmap_zero,sigmap_d,sigmap_c,kappap,pip,f_costs,kappab,Tax1,gammap),guess,options);
    
    [resid1, rl_rd_spreadp, r_l_subsidyp] = function_r_l_risky(x,betap,xip,deltap,alphap,sigma,sigma_lower_bar,taup,sigmap_zero,sigmap_d,sigmap_c,kappap,pip,f_costs,kappab,Tax1,gammap);
    
    r_l=x(1);
    
    
    
    
    % r_l = fsolve(@(r_l) -r_e ...
    %             + ( sigma_lower_bar*taup/sqrt(2*pip)...
    %             *exp(-(  ( (r_l-rl_rd_spreadp)*(1-gammap)-r_l+xip*sigma_lower_bar )/...
    %                       (sigma_lower_bar*sqrt(2)*taup)...
    %                    )^2 ...
    %                  )...
    %              + 0.5*( r_l-sigma_lower_bar*xip-(1-gammap)*(r_l-rl_rd_spreadp) )...
    %                   *(1-erf( ( (r_l-rl_rd_spreadp)*(1-gammap)-r_l+xip*sigma_lower_bar )/...
    %                       (sigma_lower_bar*sqrt(2)*taup)...
    %                          )...
    %                     )...
    %                 )/gammap...
    %         ,1/betap);
    
    r_lr = r_l*(1-r_l_subsidyp) - xip;
    
    % by parametric restriction
    r_d = r_l-rl_rd_spreadp;
    
    %
    l_upper_bar = ((r_l*(1-r_l_subsidyp)-(1-deltap))/alphap)^(1/(alphap-1));
    
    e_upper_bar = gammap*l_upper_bar;
    
    d_upper_bar = l_upper_bar-e_upper_bar;
    
    
    
    
    %k_r = sigma_lower_bar*l;
    k_r = sigma*l_upper_bar;
    
    
    k = (1-sigma)*l_upper_bar;
    
    k_ss = k;
    
    invest = deltap*(k+k_r);
    
    h = k/(k+k_r);
    
    h_r = 1 - h;
    
    y = k^(alphap)*h^(1-alphap);
    
    y_r = k_r^alphap*h_r^(1-alphap) + 0.5*k_r - 0.5*(1 + 2*xip)*k_r;
    
    ytot = y+y_r;
    
    Tax = -( -sigma*l_upper_bar*taup/sqrt(2*pip) ...
        *exp(-(  ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
        (sigma*sqrt(2)*taup)...
        )^2 ...
        ) ...
        + 0.5*( r_l*l_upper_bar-r_d*d_upper_bar-sigma*l_upper_bar*xip-f_costs*l_upper_bar ) ...
        *(1+erf( ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
        (sigma*sqrt(2)*taup)...
        ) ...
        ) ...
        )/(y+y_r);
end

w = (1-alphap)*y/h;

a = 1;

dummy_instr = 1;


%lambda_c = (sigmap_zero*d^(-sigmap_d))/(1-r_d*betap);

%c = (lambda_c/(1-betap*kappap))^(-1/sigmap_c)/(1-kappap);

% Tax = (ytot - c - invest)/ytot;

c = y+y_r-invest-e_upper_bar*kappab;

c_ss = c;

lambda_c = (1-betap*kappap)*(c-kappap*c)^(-sigmap_c) ;

r_rf = 1/betap;

r_l_lead = r_l;

% r_e_lead = r_e;

% from 25, need to pin down sigma_0 to back up the choice of spread.
%sigmap_zero=( lambda_c*(1-betap*r_d) )/d^(-sigmap_d);

r_d_ss = r_d;

r_l_ss = r_l;

d_ss = 0;

l_ss = 0;

ytot_ss = ytot;

k_r_lag = k_r;

sigma_lag = sigma;

dp = r_d*d;


tot_ret = r_l*k+r_lr*k_r;

ktot = k+k_r;

ktot_ss = ktot;


Tax_lead = Tax;

profits = ( sigma*l_upper_bar*taup/sqrt(2*pip) ...
    *exp(-(  ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    )^2 ...
    ) ...
    + 0.5*( r_l*l_upper_bar-r_d*d_upper_bar-sigma*l_upper_bar*xip-f_costs*l_upper_bar ) ...
    *(1-erf( ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    ) ...
    ) ...
    )-(1+kappab)*e_upper_bar;

profits_ss = profits;

expected_dividend_upper_bar = ( sigma*l_upper_bar*taup/sqrt(2*pip) ...
    *exp(-(  ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    )^2 ...
    ) ...
    + 0.5*( r_l*l_upper_bar-r_d*d_upper_bar-sigma*l_upper_bar*xip-f_costs*l_upper_bar ) ...
    *(1-erf( ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    ) ...
    ) ...
    );

expected_dividend_ss = 0;

expected_spread = betap*0.5*(r_l-(1-gammap)*r_d-f_costs)*(1-erf( ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    ) ...
    )*l_upper_bar;
gambling_bonus = betap*sigma*l_upper_bar*taup/sqrt(2*pip) ...
    *exp(-(  ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    )^2 ...
    );

gambling_loss = betap*0.5*sigma*l_upper_bar*xip*(1-erf( ( f_costs+r_d*(1-gammap)-r_l+xip*sigma )/...
    (sigma*sqrt(2)*taup)...
    ) ...
    );  


q = 1;

ginvest = invest;
ninvest = invest;

gamma_var = gammap;

utilc = ((c-kappap*c)^(1-sigmap_c)-1)/(1-sigmap_c);

welfc = utilc/(1-betap);

utild = sigmap_zero*(d_upper_bar^(1-sigmap_d)-1)/(1-sigmap_d);

welfd = utild/(1-betap);

Util1 = utilc+utild;

Util2 = Util1;

welf = Util1/(1-betap);

epsilon_star = (f_costs+r_d*(1-gamma_var)-r_l)/sigma_bar;

epsilon_star_expected = epsilon_star;

non_defaulted = 0.5 - 0.5*erf((epsilon_star+xip)/(taup*sqrt(2)));

non_defaulted_expected = non_defaulted;

return_non_defaulted = r_l - r_d*(1-gamma_var) - xip*sigma_bar - f_costs;

return_non_defaulted_expected = return_non_defaulted;

omega_1 = return_non_defaulted*non_defaulted*l_upper_bar;

omega_1_expected = omega_1;

omega_2 = sigma_bar*taup/sqrt(2*pip)*l_upper_bar*exp(-(((epsilon_star+xip)/(taup*sqrt(2)))^2));

omega_2_expected = omega_2;

v = profits_ss/(1-betap);

ds = d_upper_bar;

l = 0;

% combine 1) and 2) below then solve for e_upper_bar and l_upper_bar
% notice that, from 2) gamma_var = e_upper_bar/l_upper_bar


%1)
% (1+kappab)*gamma_var = betap*lambda_c(1)/lambda_c*( sigma_upper_bar/q*taup/sqrt(2*pip)
%             *exp(-(  ( F_costs/l*q+f_costs*q+r_d*(1-gamma_var)*q-r_l(1)*q+xip*sigma_upper_bar)/
%                         (sigma_upper_bar*sqrt(2)*taup) 
%                    )^2 
%                  )
%              + 0.5*( r_l(1)-sigma_upper_bar*xip/q-(1-gamma_var)*r_d-f_costs )
%                   *(1-erf( ( F_costs/l*q+f_costs*q+r_d*(1-gamma_var)*q-r_l(1)*q+xip*sigma_upper_bar )/
%                         (sigma_upper_bar*sqrt(2)*taup) 
%                          )  
%                     )
%            );


%2)
e = 0;

d = 0;

expected_dividend = 0;

%chi2_lower_bar, chi2_upper_bar, mu_lower_bar, mu_upper_bar, r_e, r_e_upper_bar

chi2_upper_bar = 0;

%chi2_upper_bar = 0;

chi2_lower_bar = (1+kappab)*gammap - betap*( sigma_lower_bar/q*taup/sqrt(2*pip)...
            *exp(-(  ( f_costs*q+r_d*(1-gammap)*q-r_l*q+xip*sigma_lower_bar)/...
                        (sigma_lower_bar*sqrt(2)*taup)... 
                   )^2 ... 
                 )...
             + 0.5*( r_l-sigma_lower_bar*xip/q-(1-gammap)*r_d-f_costs )...
                  *(1-erf( ( f_costs*q+r_d*(1-gammap)*q-r_l*q+xip*sigma_lower_bar )/...
                        (sigma_lower_bar*sqrt(2)*taup)... 
                         )...  
                    )...
           );

%r_e_upper_bar = r_e;       

r_e = ( sigma_lower_bar/q*taup/sqrt(2*pip)...
            *exp(-(  ( f_costs*q+r_d*(1-gammap)*q-r_l*q+xip*sigma_lower_bar)/...
                        (sigma_lower_bar*sqrt(2)*taup) ...
                   )^2 ...
                 )...
             + 0.5*( r_l-sigma_lower_bar*xip-(1-gammap)*r_d-f_costs )...
                  *(1-erf( ( f_costs*q+r_d*(1-gammap)*q-r_l*q+xip*sigma_lower_bar )/...
                        (sigma_lower_bar*sqrt(2)*taup) ...
                         )  ...
                    )...
           )/gammap;

r_e_lead = r_e;
       
mu_upper_bar = 0;

%mu_upper_bar = 0;

mu_lower_bar = (1+kappab)*lambda_c-betap*lambda_c*r_e;

r_e_lower_bar_lead = r_e;

r_e_upper_bar_lead = r_e_upper_bar;

mu = 1;
a_lag = a;
tau_var = taup;
kap_quality = 1;

chi2_upper_bar_ss = chi2_upper_bar;

chi2_upper_bar_proposed = chi2_upper_bar;
chi2_lower_bar_proposed = chi2_lower_bar;

isp = 1;

ltot = l_upper_bar+l;
ltot2ytot = ltot/ytot;

ltot2ytot_ss = ltot2ytot;

banking_spread = r_l_lead- r_d*(1-gamma_var); 
banking_spread_ss = banking_spread;
% solve this for gamma
% expected_dividend_upper_bar = (betap)*( sigma_lower_bar*1/gamma_var/q*taup/sqrt(2*pip)...
%             *exp(-(  ( f_costs*q+r_d*(1-gamma_var)*q-r_l*q+xip*sigma_lower_bar )/...
%                         (sigma_lower_bar*sqrt(2)*taup)...
%                    )^2 ...
%                  )...
%              + 0.5*( r_l/gamma_var-r_d*(1/gamma_var-1)-sigma_lower_bar/gamma_var*xip/q-f_costs/gamma_var)...
%                   *(1-erf( ( f_costs*q+r_d*(1-gamma_var)*q-r_l*q+xip*sigma_lower_bar )/...
%                         (sigma_lower_bar*sqrt(2)*taup)...
%                          )...
%                     )...
%                                       )-(1+kappab);
                                  


