function [ys,params,check] = Bench_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = Argentina_steadystate(ys,exo,M_,options_)
% computes the steady state for Argentina.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
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)


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


%% Enter model equations here
options=optimset(); % set options for numerical solver

    % DISTURBANCES AND AUXILIARY VARIABLES
    prefshock = 0;
    UIP = 0;
    premA = 0;
    ssmus = 1 - (1 - ssmu) * M / Ms;
    mu1 = ssmu;
    mu2 = 1 - ssmu;
    mus2 = 1 - ssmus;
    mus1 = 1 - mus2;
    size = M / (M + Ms);
    sizes = Ms / (M + Ms);
    pmark = sig / (sig - 1);
    pmarks = pmark;
    rmark  = sigR/(sigR - 1);
    rmarks = sigR/(sigR - 1);
    wmark = phi / (phi - 1);
    wmarks = phi / (phi - 1);
    z = 0;
    zs = 0;

    % PRICES AND INFLATION
    rer   = 0;
    delS  = 0;
    tot   = 0;
    cpin = 0;
    cpins = 0;
    mcr  = -log(rmark);
    mcrs = -log(rmarks);
    infl   = 0;
    infls  = 0;
    infh  = 0;
    inffs = 0;
    winfl  = 0;
    winfls = 0;
    wn  = 0;
    wns = 0;
    DPh   = 1;
    DPfs  = 1;

    % NUMERICAL STARTING VALUES
    ph = -0.205322;
    pfs = -0.181061;

    % STEADY STATE REDUCED FORM EQUATIONS
    int  = (1/betta) - 1;
    ints = (1/bettas) - 1;
    r   = log( delta + int );
    rs  = log( delta + ints );
    % IntPr = int - ints;
    w = log(1 - nu) + (1/(1-nu)) * ( ph + nu * log(nu) - nu * r - log(pmark) );
    ws = log(1 - nu) + (1/(1-nu)) * ( pfs + nu * log(nu) - nu * rs - log(pmarks) );
    mc = nu * r + (1-nu) * w - nu * log(nu) - (1-nu) * log(1-nu);
    mcs = nu * rs + (1-nu) * ws - nu * log(nu) - (1-nu) * log(1-nu);
    klr  = w - r + log(nu) - log(1 - nu);
    klrs = ws - rs + log(nu) - log(1 - nu);
    ilr  = log(delta) + klr;
    ilrs = log(delta) + klrs;
    ylr  = nu * klr + log(M) + z;
    ylrs = nu * klrs + log(Ms) + zs;
    clr  = log( (1/M) * ( exp(ylr)*(1 - ssmu*gshare) - exp(ylrs)*(1-ssmus)*gshares ) - exp(ilr) );
    clrs = log( (1/Ms) * ( exp(ylrs)*(1 - ssmus*gshares) - exp(ylr)*(1-ssmu)*gshare ) - exp(ilrs) );
    l  = (w  - log(pmark) - theta  * clr) /(theta  + chi);
    ls = (ws - log(pmark) - thetas * clrs)/(thetas + chi);
    k = klr + l;
    ks = klrs + ls;
    invt = ilr + l;
    invts = ilrs + ls;
    invh = log(mu1) + invt;
    invf = log(mu2) + invt;
    invfs = log(mus1) + invts;
    invhs = log(mus2) + invts;
    y = ylr + l;
    ys = ylrs + ls;
    yh = log(mu1) + y;
    yf = log(mu2) + y;
    yfs = log(mus1) + ys;
    yhs = log(mus2) + ys;
    gdp = log( exp(ph-mcr+y) + (1-exp(-mcr))*exp(ph+yhs) + (exp(-mcr)-1)*exp(rer+pfs+yf) );
    gdps = (pfs - mcrs) + ys;
    c = clr + l;
    cs = clrs + ls;
    ch = log(mu1) + c;
    cf = log(mu2) + c;
    cfs = log(mus1) + cs;
    chs = log(mus2) + cs;
    lam  = -theta * c;
    lams = -thetas * cs;
    ksi  = exp(lam);
    ksis = exp(lams);
    g = log(gshare) + gdp - log(M);
    gs = log(gshares) + gdps - log(Ms);
    gh = log(mu1) + g;
    gf = log(mu2) + g;
    gfs = log(mus1) + gs;
    ghs = log(mus2) + gs;

    wb   = exp(l * (1 + chi)) / (1 - omega * betta);
    wbs  = exp(ls * (1 + chi)) / (1 - omegas * betta);
    wa   = exp(lam + w + l) / (1 - omega * betta);
    was  = exp(lams + ws + ls) / (1 - omegas * betta);
    pbh  = exp(lam + sig*ph + mc + yh) / (1 - alphap * betta);
    pah  = exp(lam + sig*ph + yh) / (1 - alphap * betta);
    pbfs = exp(lams + sig*pfs + mcs + yfs) / (1 - alphaps * betta);
    pafs = exp(lams + sig*pfs + yfs) / (1 - alphaps * betta);
    phn  = ph;
    pfsn = pfs;

    DW  = M;
    DWs = Ms;
    AL  = DW  * exp(l  * (1 + chi));
    ALs = DWs * exp(ls * (1 + chi));
    u  = M  * c  - AL  / (1 + chi);
    us = Ms * cs - ALs / (1 + chi);
    v  = u  / (1 - betta);
    vs = us / (1 - bettas);
    disl  = -AL  / ((1 + chi) * (1 - betta));
    disls = -ALs / ((1 + chi) * (1 - bettas));

    wrate = exp(w - ws - rer);
    crate = exp(c - cs - rer);
    relc  = c - cs;
    rely  = y - ys;
    relz = z - zs;
    rell = l - ls;
    gr  = M * exp(g - y);
    grs = Ms * exp(gs - ys);
    cr  = M * exp(c - y);
    crs = Ms * exp(cs - ys);
    ir  = M * exp(invt - y);
    irs = Ms * exp(invts - ys);
    kr  = M * exp(k - y)/4;
    krs = Ms * exp(ks - ys)/4;


    yret  = log(exp(c)  + exp(invt)  + exp(g));
    yrets = log(exp(cs) + exp(invts) + exp(gs));
    pbc = exp( lam + mcr + yret)/ (1 - alphar * betta);
    pac = exp( lam + yret ) / (1 - alphar * betta);
    pbcs = exp( lams + mcrs + yrets) / (1 - alphars * bettas);
    pacs = exp( lams + yrets) / (1 - alphars * bettas);
    
    tr  = -exp(g);
    trs = -exp(gs);
    trr  = -M * tr / exp(y);
    trrs = -Ms * trs / exp(ys);

    nfwr = - nfdr;
    nfw = nfwr*(4*exp(gdp));
    nfd = - nfw;
    dr = nfd / (4*exp(gdp));
    IntPr = (1/betta - 1/bettas)+ exp( risk * (nfd/(4*exp(gdp)) - nfdr) ) - 1 + UIP;
    ca =    nfw * (ints + IntPr) * exp(delS) / exp(infl) + Ms * exp(ph) * (exp(chs) + exp(invhs) + exp(ghs)) - ( exp(pfs + rer) * M * (exp(cf) + exp(invf) + exp(gf)) );
    ex = Ms * exp(ph) * (exp(chs) + exp(invhs) + exp(ghs));
    im = ( exp(pfs + rer) * M * (exp(cf) + exp(invf) + exp(gf)) );
    tbal =  ex - im;
    car = ca / exp(y);
    expend = log(M*exp(c) + M*exp(invt) + M*exp(g) + ex - im);
 
    ikr  = exp(invt - k) - delta;
    ikrs = exp(invts - ks) - delta;
    err = exp(c) + exp(invt) - tr + nfw - (exp(gdp) + (1 + ints + IntPr) * nfw * exp(delS) / exp(infl));
    
%% end own model equations

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

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