function [ys,params,check] = Model_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% 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
%   - 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)

% Copyright (C) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

%% Step 0: initialize indicator and set options for numerical solver
check=0;
options= optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

%% Step 1: read out parameters to access them with their name

for ii = 1:M_.param_nbr                                 
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);         %    Get the value of parameter i.
end                                                                     % End of the loop.  




%% Step 2: Enter model equations here


% values of the parameters I use here: 
b_eta=0.995;        % discount factor
s_igma=2;          % inverse intertemporal elasticities of substitution
eta=2;             % inverse labor supply elasticity
g_amma=0.23;       % long run share of foreign consumption
a=5;               % elasticity of domestic consumption
theta=6;           % implies a mark up of 1.2
omega=0.75;        % domestic Calvo parameter
rho_y=0.9;         % persistence foreign demand shock
rho_d=0.9;         % persistence terms of trade shock
sigma_y=0.01;
sigma_d=0.01;

% log-lin variables in SS:

ex_h=0;
y_h=0;
im_h=0;
pi_c_h=0;
a1_h=0;
a2_h=0;
temp_h=0;
ystar_h=0;
delta_h=0; 
y_obs =0;
pi_c_obs = 0;
im_obs = 0; 
ex_obs = 0; 
delta_obs = 0; 
ystar_obs = 0; 
%nex_obs=0;
%ey=0; % Ireland 

% variables values in SS: 
delta=1;
temp=1;
pi_c=1;
pi_h=1;
psi=(theta-1)/theta;
y=psi^(1/(eta+s_igma));
w=psi;
ystar=y;
c=y;
n=y;
ex=g_amma*ystar;
im=g_amma*ystar;
a1=y*psi/(1-omega*b_eta);
a2=y/(1-omega*b_eta);




a51=-(a1*b_eta*g_amma*omega*temp*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1) - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(theta + 1)*(delta^(1 - a)*g_amma - g_amma + 1))/(temp^a*(g_amma*temp^(1 - a) - g_amma + 1)^2);
a52=(a1*b_eta*delta*omega*s_igma*ystar*ystar^(s_igma - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta + 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1)*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma));
a53=delta*((a1*b_eta*omega*(ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta + 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma)) + (a1*b_eta*g_amma*omega*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1) - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(theta + 1))/(delta^a*(g_amma*temp^(1 - a) - g_amma + 1)));
a57=pi_c*((a1*b_eta*omega*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta + 1))/pi_c^2 - (a1*b_eta*omega*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1))*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(theta + 1))/pi_c);
a58=-(a1*b_eta*omega*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta + 1))/pi_c;
a61=-(a2*b_eta*g_amma*omega*temp*theta*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1) - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta - 1)*(delta^(1 - a)*g_amma - g_amma + 1))/(temp^a*(g_amma*temp^(1 - a) - g_amma + 1)^2);
a62=(a2*b_eta*delta*omega*s_igma*ystar*ystar^(s_igma - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1)*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma));
a63=delta*((a2*b_eta*omega*(ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma)) + (a2*b_eta*g_amma*omega*theta*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1) - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta - 1))/(delta^a*(g_amma*temp^(1 - a) - g_amma + 1)));
a67=pi_c*((a2*b_eta*omega*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta)/pi_c^2 - (a2*b_eta*omega*theta*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1))*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta - 1))/pi_c);
a69=-(a2*b_eta*omega*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta)/pi_c;
a71=temp;
a82=ystar;
a93=delta;

b12=ystar*(delta^a*g_amma - (delta*ystar^(s_igma - 1)*(g_amma - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1)*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))/(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1)));
b13=-delta*(((ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*(g_amma - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1))/(s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1))) - a*delta^(a - 1)*g_amma*ystar + (a*g_amma*(g_amma - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))/(delta^a*(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1) + 1)));
b14=-y;
b22=delta^a*g_amma*ystar;
b23=a*delta*delta^(a - 1)*g_amma*ystar;
b25=-ex;
b32=(delta*delta^(1 - a)*g_amma*ystar*ystar^(s_igma - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1)*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))/(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1));
b33=delta*((delta^(1 - a)*g_amma*(ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1))/(s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1))) - (g_amma*(a - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))/(delta^a*(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1))) + (a*delta^(1 - a)*g_amma^2*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))/(delta^a*(delta^(1 - a)*g_amma - g_amma + 1)^(a/(a - 1) + 1)));
b36=-im;
b41=(g_amma*omega*pi_c*temp*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1) - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta - 2)*(theta - 1)*(delta^(1 - a)*g_amma - g_amma + 1))/(temp^a*(g_amma*temp^(1 - a) - g_amma + 1)^2);
b43=-(delta*g_amma*omega*pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1) - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta - 2)*(theta - 1))/(delta^a*(g_amma*temp^(1 - a) - g_amma + 1));
b47=omega*pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1))*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta - 2)*(theta - 1);
b48=(a1*theta*(omega - 1))/(a2*((a1*theta)/(a2*(theta - 1)))^theta);
b49=-(a1*theta*(omega - 1))/(a2*((a1*theta)/(a2*(theta - 1)))^theta);
b52=ystar*(delta*s_igma*y*y^eta*ystar^(s_igma - 1)*((delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))^(s_igma - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1) + (a1*b_eta*delta*omega*s_igma*ystar^(s_igma - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta + 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma)*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma + 1)));
b53=delta*((g_amma*y*y^eta*((delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))^s_igma)/(delta^a*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) + 1)) + (y*y^eta*(ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*((delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))^(s_igma - 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma - 1))/(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) + (a1*b_eta*omega*(ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^(theta + 1)*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma + 1)));
b54=y*((y^eta*((delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))^s_igma)/(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) + (eta*y*y^(eta - 1)*((delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))^s_igma)/(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)));
b58=-a1;
b62=(a2*b_eta*delta*omega*s_igma*ystar*ystar^(s_igma - 1)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma)*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma + 1));
b63=(a2*b_eta*delta*omega*(ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)) - (delta*g_amma*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1) - 1))/delta^a)*(pi_c*((delta^(1 - a)*g_amma - g_amma + 1)/(g_amma*temp^(1 - a) - g_amma + 1))^(1/(a - 1)))^theta*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma))/(pi_c*(delta*ystar^s_igma*(delta^(1 - a)*g_amma - g_amma + 1)^(1/(a - 1)))^(1/s_igma + 1));
b64=y;
b69=-a2;
b73=delta;
b82=rho_y*ystar;
b93=delta*delta^(rho_d - 1)*rho_d;





%% Step 3: Update parameters and variables

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

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


