% Staeady State Model for Fiscak Multiplier Model

function [ys,check] = keynes_fiscal_mult_steadystate(ys,exo)
% function [ys,check] = NK_baseline_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 of 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;

sub1=0;
tfp=1;
nnu=1;
tc = tcbar;
tw = twbar;
tk=tkbar;
tz = tzbar;
pz=1;
nxy=nxybar;
ig=igbar;
cg=cgbar;
pii=piibar;
slr= slrbar;
crr=crrbar;
RG=piibar/betta;
RD = pii/betta;
RE = (RD - slr*RG - crr)/(1-slr-crr);
q = 1+ ppsi*((RE/RD)-1);
piistar = ((pii^(1-xii_y) - gammap*pii^(iotap*(1-xii_y)))/(1-gammap))^(1/(1-xii_y));
pstar = piistar/pii;
pf = ((xii_y-1)/xii_y)*pstar;
Vp = ((1-gammap)*(piistar/pii)^(-xii_y))/(1-gammap*pii^((1-iotap)*xii_y));
s = (cchi + (1-cchi)*(1+tz)^(1-zzeta))^(1/(1-zzeta));
rk = (1/betta - 1 + ddelta)*(s/(1-tk));
kg = ig/ddelta;
ihg = cchi*ig*s^zzeta;
izg = (1-cchi)*ig*(s/(1+tz))^zzeta;
mpk = q*rk;
kfp_v = (aalpha_1*(1+tz))/(aalpha_2*rk);
l_v = ((mpk*kfp_v^(1-aalpha_1))/(aalpha_1*pf*tfp*kg^aalpha_g))^(1/(1-aalpha_1-aalpha_2));
mpl = (1-aalpha_1-aalpha_2)*tfp*pf*kg^aalpha_g*kfp_v^aalpha_1* l_v^(-aalpha_1-aalpha_2); 
w = mpl/q;
wstar = w*((1-gammaw*pii^((iotaw-1)*(1-xii_l)))/(1-gammaw))^(1/(1-xii_l));
wh = ((xii_l-1)/xii_l)*wstar;
Vw = ((1-gammaw)*(wstar/w)^(-xii_l))/ (1-gammaw*pii^((1-iotaw)*xii_l));

%  Labor Supply of Non-Ricardian households

xx0=[0.21 , 0.04];

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

xx= fsolve(@consumption, xx0,  options, wh,  kg, l_v, kfp_v, s,  ig, Vp, Vw, mmunr, kkappa, tcbar, eepsilon, oomega, cgbar, mmur,   twbar,   ttheta, ddelta, aalpha_g, aalpha_1, aalpha_2, nxybar);
cpr=xx(1);
cpnr=xx(2);

cnrhat = cpnr^((eepsilon-1)/eepsilon)+ oomega*cgbar^((eepsilon-1)/eepsilon);
lambdanr= (cpnr^(-1/eepsilon))/((1+tc)*cnrhat);
lhnr = 1-(mmunr/(kkappa*wh*lambdanr));


% Loop for the Labor Supply of Ricardian households

    
crhat = cpr^((eepsilon-1)/eepsilon)+ oomega*cg^((eepsilon-1)/eepsilon);
lambdar = (nnu/((1+tc)*crhat))*cpr^(-1/eepsilon);
lhr = 1- (mmur/(wh*(1-tw)*lambdar));
lh = ttheta*lhr + (1-ttheta)*kkappa*lhnr;
l= lh/Vw;
v= (1/l_v)*l;
tilde_kp = kfp_v*v; % Firms Demand for Capital
kp = tilde_kp/ttheta; % Capital Stock per Ricardian household
ip= ddelta*kp;% Investment per Ricardian household
yf = tfp*kg^aalpha_g*tilde_kp^ aalpha_1*v^ aalpha_2*l^( 1-aalpha_1- aalpha_2);
y= yf/Vp;
nx= nxy*y;

be= (ppsi*(w*l+rk*tilde_kp+v*(1+tz)))/pii;
d = (be/((1-slr-crr)*ttheta));
bc = crr*ttheta*d;
bg= slr*ttheta*d;
ihp = cchi*ip*s^zzeta;
izp= (1-cchi)*ip*(s/(1+tz))^zzeta;
ystar= y*pstar^(-xii_y);
lstar = l*wstar^(-xii_l);
cp = ttheta*cpr+(1-ttheta)*cpnr;
cr = crhat^(eepsilon/(eepsilon-1));
cnr = cnrhat^(eepsilon/(eepsilon-1));
f1= (lambdar*ystar*pstar)/(1-gammap*betta);
f2= (lambdar*ystar*pf)/(1-gammap*betta);
g1= (lambdar*lstar*wstar)/(1-gammaw*betta);
g2= (lambdar*lstar*wh)/(1-gammaw*betta);
ln_y = log(y);
 ln_cp = log(cp);
ln_cg =log(cg);
 ln_ig= log(ig);
ln_ip= log(ip);



for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

M_.params(28) = y;
M_.params(29) = bg;



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
end

