% Staeady State Model for Fiscak Multiplier Model

function [ys,check] = keynes_multiplier_steadystate(ys,exo)
% function [ys,check] = keynes_multiplier_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 if not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator

check = 0;

tfps=1;
nnus=1;
tcs = tcbar;
tws = twbar;
tks=tkbar;
tzs = tzbar;
pzs=1;
nxys=nxybar;
igs=igbar;
cgs=cgbar;
piis=piibar;
slrs= slrbar;
crrs=crrbar;
RGs=piibar/betta;
RDs = piis/betta;
rds=1/betta;
REs = (RDs - slrs*RGs - crrs)/(1-slrs-crrs);
qs = 1+ ppsi*((REs/RDs)-1);
piistars = ((piis^(1-xii_y) - gammap*piis^(iotap*(1-xii_y)))/(1-gammap))^(1/(1-xii_y));
pstars = piistars/piis;
pfs = ((xii_y-1)/xii_y)*pstars;
Vps = ((1-gammap)*(piistars/piis)^(-xii_y))/(1-gammap*piis^((1-iotap)*xii_y));
ss = (cchi + (1-cchi)*(1+tzs)^(1-zzeta))^(1/(1-zzeta));
rks = (1/betta - 1 + ddelta)*(ss/(1-tks));
kgs = igs/ddelta;
ihgs = cchi*igs*ss^zzeta;
izgs = (1-cchi)*igs*(ss/(1+tzs))^zzeta;
mpks = qs*rks;
kfp_v = (aalpha_1*(1+tzs))/(aalpha_2*rks);
l_v = ((mpks*kfp_v^(1-aalpha_1))/(aalpha_1*pfs*tfps*kgs^aalpha_g))^(1/(1-aalpha_1-aalpha_2));
mpls = (1-aalpha_1-aalpha_2)*tfps*pfs*kgs^aalpha_g*kfp_v^aalpha_1* l_v^(-aalpha_1-aalpha_2); 
ws = mpls/qs;
wstars = ws*((1-gammaw*piis^((iotaw-1)*(1-xii_l)))/(1-gammaw))^(1/(1-xii_l));
whs = ((xii_l-1)/xii_l)*wstars;
Vws = ((1-gammaw)*(wstars/ws)^(-xii_l))/ (1-gammaw*piis^((1-iotaw)*xii_l));

xx0=[.21 , .04];

options = optimoptions('fsolve', 'display', 'off');

[x, fval, exitflag]= fsolve(@consumption_bay, xx0,  options, whs,  kgs, l_v, kfp_v, ss,  igs, Vps, Vws, mmunr, kkappa, tcbar, eepsilon, oomega, cgbar, mmur,   twbar,   ttheta, ddelta, aalpha_g, aalpha_1, aalpha_2, nxybar);

if exitflag <1 
    %indicate the SS computation was not sucessful; 
     check=1; %set failure indicator return; %return without updating steady states
     return;
end 
    
    
if  isreal(x(1)) ==0 
    check=1; %complex value; 
    return; %return without updating steady states 
end 

if  x(1) <=0 
    check=1; %negative value; 
    return; %return without updating steady states 
end 

if  isreal(x(2)) ==0 
    check=1; %complex value; 
    return; %return without updating steady states 
end

if  x(2) <=0 
    check=1; %negative value; 
    return; %return without updating steady states 
end 


cprs=x(1);
cpnrs=x(2);

cnrhats = cpnrs^((eepsilon-1)/eepsilon)+ oomega*cgs^((eepsilon-1)/eepsilon);
lambdanrs= (cpnrs^(-1/eepsilon))/((1+tcs)*cnrhats);
lhnrs = 1-(mmunr/(kkappa*whs*lambdanrs));

if  lhnrs <=0 
    check=1;  %return without updating steady states 
    return;
end

if  lhnrs >=1  
    check=1;  
    return; %return without updating steady states 
end


crhats = cprs^((eepsilon-1)/eepsilon)+ oomega*cgs^((eepsilon-1)/eepsilon);
lambdars = (nnus/((1+tcs)*crhats))*cprs^(-1/eepsilon);
lhrs = 1- (mmur/(whs*(1-tws)*lambdars));

if  lhrs <=0 
    check=1;  %return without updating steady states 
    return;
end

if  lhrs >=1  
    check=1; 
    return; %return without updating steady states 
end



lhs = ttheta*lhrs + (1-ttheta)*kkappa*lhnrs;
ls= lhs/Vws;
vs= (1/l_v)*ls;
tilde_kps = kfp_v*vs; % Firms Demand for Capital
kps = tilde_kps/ttheta; % Capital Stock per Ricardian household
ips= ddelta*kps;% Investment per Ricardian household
yfs = tfps*kgs^aalpha_g*tilde_kps^ aalpha_1*vs^ aalpha_2*ls^( 1-aalpha_1- aalpha_2);
ys1= yfs/Vps;
nxs= nxys*ys1;

bes= (ppsi*(ws*ls+rks*tilde_kps+vs*(1+tzs)))/piis;
ds = (bes/((1-slrs-crrs)*ttheta));
bcs = crrs*ttheta*ds;
bgs= slrs*ttheta*ds;
ihps = cchi*ips*ss^zzeta;
izps= (1-cchi)*ips*(ss/(1+tzs))^zzeta;
ystars= ys1*pstars^(-xii_y);
lstars = ls*wstars^(-xii_l);
cps = ttheta*cprs+(1-ttheta)*cpnrs;
crs = crhats^(eepsilon/(eepsilon-1));
cnrs = cnrhats^(eepsilon/(eepsilon-1));
f1s= (lambdars*ystars*pstars)/(1-gammap*betta);
f2s= (lambdars*ystars*pfs)/(1-gammap*betta);
g1s= (lambdars*lstars*wstars)/(1-gammaw*betta);
g2s= (lambdars*lstars*whs)/(1-gammaw*betta);

sales_taxs=(tcs*(ttheta*cprs + (1-ttheta)*cpnrs)); 
wage_taxs= (ttheta*tws*whs*lhrs);
capital_taxs= ( ttheta* rks*kps*tks);
custom_taxs = tzs*(vs+ttheta*izps+ izgs);

y=log(ys1);
gdp=log(ys1);
yf=log(yfs);
ystar=log(ystars);
cp =log(cps);
cr=log(crs);
cnr =log(cnrs);
cpr =log(cprs);
cpnr =log(cpnrs);  
ip =log(ips);
ihp =log(ihps);
izp =log(izps);
kp =log(kps);
tilde_kp =log(tilde_kps);
kgr=0;
kg =log(kgs);
cg =log(cgs);
ig =log(igs);
ihg =log(ihgs);
izg =log(izgs);
nx =nxs;
nxy = nxys;
lh =log(lhs);
lhr =log(lhrs);
lhnr =log(lhnrs);
l =log(ls);
lstar =log(lstars);
v =log(vs);
tfp =log(tfps);
nnu =log(nnus);
RG =log(RGs);
RD =log(RDs);
rd=log(rds);
RE =log(REs); 
pii =log(piibar);
piistar =log(piistars);
rk =log(rks);
w =log(ws);
wh =log(whs);
wstar =log(wstars);
pf =log(pfs);
pstar =log(pstars);
pz =log(pzs);
q =log(qs);
Vp =log(Vps);
Vw =log(Vws);
s = log(ss);
be =log(bes);
bc =log(bcs);
bg =log(bgs);
d =log(ds);
crr =log(crrs);
slr =log(slrs);
tc =log(tcs);
tw =log(tws);
tk =log(tks);
tz =log(tzs);
lambdar =log(lambdars);
lambdanr =log(lambdanrs);
f1 =log(f1s); 
f2 =log(f2s);
g1 =log(g1s);
g2=log(g2s);
crhat=log(crhats);
cnrhat =log(cnrhats);

% Relation to observables

sales_tax = log(sales_taxs);
capital_tax = log(capital_taxs); 
wage_tax = log(wage_taxs);
custom_tax = log(custom_taxs);



NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end


M_.params(29)=ys1; 
M_.params(30)=bgs; 

end

