function [ys,check] = OBC_DSCES_392016_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;
%%
%%Model
load par.mat
%%Load some variables%%%%%%%%%%%%%%%%%%%
V=theta*u;%Vacancies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms MC b chi
w=k/(B*theta^(-lambda))*(1/0.4);
G(1)=MC-(1/Z)*(w+k*((1-betta*(1-delta)*(1-s))/(B*theta^(-lambda))));
G(2)=w-(chi*(MC*Z+betta*k*theta)+(1-chi)*b);
G(3)=b-0.95*w;
[MC b chi]=vpasolve(G(1:3),[MC b chi])
p=MC*1/rho;
MC=double(MC);b=double(b);chi=double(chi);p=double(p);
%%%%%%%%%%%%%%
%%Define structure of parameters
param.MC=MC;
param.b=b;
param.chi=chi;
param.w=w;
param.p=p;
param.eps=0.1;
param.rho=rho;
param.betta=betta;
param.delta=delta;
param.u=u;
param.s=s;
param.theta=theta;
param.V=V;
param.lambda=lambda;
param.C=B;
param.Z=Z;
param.Psi=Psi;
param.r=r;
param.nu=nu;
param.R=R;
param.omega=omega;
param.k=k;
param.A=A;
if R > Rbar % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form
    check=1; %set failure indicator
    error('Error. R must be less than Rbar')
end

%%
%%%%%%%Calibration Targets%%%%%%%%%%%%%%%%
%Average hiring cost is 40% of wages
%rho chosen to satisfy 30% markup
%b=0.95 w (Hagedorn Manovskii)
%chi implied from b=0.95w
%gam chosen to satify equilibrium u=0.051, theta=0.51
%A chosen to satisfy Ra/cons=0.816 (M2 share of consumption)
%eps?
%omega?
%k? Satisfy gam= 5.94*operating cost
    %operating cost=nbar*(w+k*s)%Barseghyan and DiceCio

   

%w,MC,chi,b,Q,Qt,qm,qc,qt,a,at,d,S,nbar,Y,cons
%p,A,gam
init=[0.4996, 0.4010, 5.9,6.1685 ,5.4771, 0.0465, 0.4317,...
    0.2193, 0.1545, 6.1417, 1.2426, 1.2402,12.7973,rand(1),rand(1),0.24,0.2]
%load init
init=init+rand(size(init))
options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective',...
    'MaxFunEvals',9000,'Display','iter','TolFun',1*10^(-20),'TolX',1*10^(-14),...
    'MaxIter',2000)
%options=optimoptions('fsolve','Display','iter','MaxFunEvals',2000)
lb=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
ub=[Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,...
    Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]
[x,resnorm,residual,exitflag,output]=lsqnonlin(@(x) roolz(x,param),init,lb,ub,options)
Q=x(1);Qt=x(2); qm=x(3);
qc=x(4); qt=x(5); a=x(6);
at=x(7); d =x(8);
 S=x(9); nbar=x(10); Y=x(11); C=x(12);
  gam=x(13);diff=x(14);diff2=x(15);diff3=x(16);
  Omega=x(17);
  %Define remaining variables
  nm=qm;
  nc=qc;
  share=S*p*(omega*qm+(1-omega)*qc)/C;
  P=p*share+1-share;%price index
  wr=w/P;%real wage
varnames={'Q', 'Qt', 'qm', 'qc', 'qt', 'a',...
    'at' ,'d', 'S', 'nbar' ,'Y' ,'C','gam','diff','diff2','diff3','Omega'};
for i=1:length(x)
myStruct.(varnames{i}) =x(i);
end
myStruct.qstar=(A*S^((1-rho-eps)/rho)/p)^(1/eps);
myStruct.b=b; mySruct.w=w; 
myStruct.MC=MC; myStruct.chi=chi;
myStruct
oper_cost=nbar*(w+k*s);
gam/oper_cost
display('share of liquid assets')
R*a/C
display('share of debt')
d/C
%F(21)=Delta-(d+R*nu*a+y-w*nbar-k*V);
exitflag
resnorm
%% end own model equations

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

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