function [ys,check] = test_v3_steadystate(ys,exo)
global M_

%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
% 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;
%%
%% END OF THE FIRST MODEL INDEPENDENT BLOCK.


%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.
%%
%% Here the user has to define the steady state.
options=optimset(); % set options for numerical solver
% options = optimoptions('fsolve','Display','iter');
% opts.Algorithm = 'levenberg-marquardt';


% stedy state value of the variables have an (ss) to them.
% steady state shock processes (domestic)
nubeta  = 1         ;            % discount factor shock at ss (domestic)
nun     = 1         ;            % labor supply shock at ss (domestic)
nuc     = 1         ;            % taste shock at ss (domestic)
nui     = 1         ;            % investment shock at steady state (domestic)
nup     = 1         ;            % productivity shock at steady state (domestic)

% steady state values of tax processes (domestic)
tauyd   = tauydss   ;            % output tax rate at steady state (domestic)
taun    = taunss    ;            % labor tax rate at steady state (domestic)
tauk    = taukss    ;            % capital tax rate at steady state (doemstic)
taud    = taudss    ;            % tax rate on doemstic consumption at ss (domestic)
taum    = taumss    ;            % tax rate on imported consumption at ss (domestic)

% steady state values of inflation processes (domestic)
pid     = pidss     ;            % domestic inflation rate in steady state (domestic) 
pi      = pidss     ;            % cpi inflation in steady state (domestic)
pidopt  = pidss     ;            % reset price inflation in steady state (domestic)
syms psigma
eq1 = psigma == (1-tauydss)^(ej)*((1-xi)*((1+pidopt)/(1+pid))^(-ej) + xi*((1+pidss)/(1+pid))^(-ej) * psigma ); 
[A,B] = equationsToMatrix([eq1], [psigma]);
XD = linsolve(A,B);
psigma = double(XD(1,1));        % price dispersion in steady state (domestic)

% real rate of return on capital in steady state (domestic)
r = (1/(1-tauk)) * ((1/betta)  - (1- delta)); 

% marginal cost in steady state (domestic)
mc = ((ej - 1)/ej);

% tempvar K/N in steady state (domestic)
kn = (mc*alppha/r)^(1/(1-alppha));

% real wage rate in steady state (domestic) ~ done really need either k or
% n for wage rate, just the ratio.
w  = (1-alppha) * mc * (kn)^(alppha); 

% foreign variables have an extra (s) to them indicating the star in these
% variables
nubetas     = 1         ;            % discount factor shock at ss (foreign)
nuns        = 1         ;            % labor supply shock at ss (foreign)
nucs        = 1         ;            % taste shock at ss (foreign)
nuis        = 1         ;            % investment shock at steady state (foreign)
nups        = 1         ;            % productivity shock at steady state (foreign)


tauyds      = tauydsss  ;            % output tax rate at steady state (foreign)
tauns       = taunsss   ;            % labor tax rate at steady state (foreign)
tauks       = tauksss   ;            % capital tax rate at steady state (foreign)
tauds       = taudsss   ;            % tax rate on doemstic consumption at ss (foreign)
taums       = taumsss   ;            % tax rate on imported consumption at ss (foreign) 


pids        = pidsss     ;            % domestic inflation rate in steady state (foreign)
pis         = pidsss     ;            % cpi inflation in steady state (foreign)
pidopts     = pidsss     ;            % reset price inflation in steady state (foreign)

% price dispersion in steady state (foreign)
syms psigmas
eq1 = psigmas == (1-tauydsss)^(ej)*((1-xi)*((1+pidopts)/(1+pids))^(-ej) + xi*((1+pidss)/(1+pids))^(-ej) * psigmas );
[A,B] = equationsToMatrix([eq1], [psigmas]);
XF = linsolve(A,B);
psigmas = double(XF(1,1));
 
% real rate of return on capital in steady state (foreign) 
rs = (1/(1-tauks)) * ((1/betta)  - (1-delta));

% marginal cost in steady state (foreign)
mcs = ((ej - 1)/ej);

% tempvar K^*/N^* in steady state (foreign)
kns = (mc*alppha/rs)^(1/(1-alppha));

% real wage rate in steady state (foreign)
ws = (1-alppha) * mc * (kns)^(alppha);

% monetary policy at steady state (cu)  
num          = 1        ;            

% nominal interest rate in ss (cu)
% when the interest rate is different the we will have interest rate as a
% result of this equality
icuss = ((1+pi)/betta) - 1;  % icuss = ((1+pid)/betta) - 1 = ((1+pids)/betta) - 1
icu = icuss;

% terms of trade between two countries at steady state (cu)
tau = tauss; 


% tempvar yd/n = f(kn) (domestic)
ydn = nup* kn^(alppha) * (psigma/((1-tauyd)^(ej)));

% tempvar yd^*/ns = f(kn) (foreign)
ydns = nups * kns^(alppha) * (psigmas/((1-tauyds)^(ej)));

% complex parameters for private consumption (domestic)
acd = ((1-om)*(1+taud)^(1-el) + om * (1+taum)^(1-el) * tau^(1-el))^(el/(1-el));
acm = ((1-om)*(1+taud)^(1-el)*tau^(-(1-el))+om*(1+taum)^(1-el))^(el/(1-el));

% complex parameters for government consumption (domestic)
agd = ((1-omg)*(1+taud)^(1-eg)+ omg*(1+taum)^(1-eg) * tau^(1-eg))^(eg/(1-eg));
agm =  ((1-omg)*(1+taud)^(1-eg)* tau^(-(1-eg)) + omg*(1+taum)^(1-eg) )^(eg/(1-eg));

% complex parameters for private consumption (foreign)
acds = ((1-oms)*(1+tauds)^(1-el) + oms*(1+taums)^(1-el) * tau^(-(1-el)))^(el/(1-el));
acms = ((1-oms)*(1+tauds)^(1-el) * tau^(1-el) + oms*(1+taums)^(1-el))^(el/(1-el));

% complex parameters for government consumption (foreign)
agds = ((1-omgs)*(1+tauds)^(1-eg) + omgs*(1+taums)^(1-eg) * tau^(-(1-eg)))^(eg/(1-eg));
agms = ((1-omgs)*(1+tauds)^(1-eg)*tau^(1-el) + omgs*(1+taums)^(1-eg))^(eg/(1-eg));

% tempvar gdn = gd/n domestic govt consumption per unit labor (domestic)
gdn = (gss * ydn*omg*agd)/(1+taud)^(eg); 
gmn = (gss* ydn*omgs*agm)/(1+taum)^(eg); % tempvar gmn = g/n


% tempvar gdns = gd*/n* domestic govt consumption per unit labor (foreign)
gdns = (gsss * ydns * omgs *agds)/(1+tauds)^(eg); 



% tempvar gm*/n*
gmns = (gss* ydn*omgs*agms)/(1+taumss)^(eg); 

% tempvar in = i/n
in = delta*kn;

% tempvar in = i*/n*
ins = delta*kns;

% type 1 - orignial version
eq = @(X) [ ydn - ( (X(1) *mun*(1-om)*acd)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taud)^(el)) ...
    + (ydn*(1-omg)*gss*agd)/((1+taud)^(eg)) + delta * kn ...
    + ((X(2)- muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns))/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)) ...
    * (zetas/zeta) * ( (X(2)*mun*om*acms)/((X(2)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)*(1+taums)^(el)) ...
    + (gsss*ydns*omg*agms)/(1+taums)^(eg) )) ,...
    ydns - ( (X(2)*mun*(1-oms)*acds)/((X(2)-muc*nucs)^(-ec)*(1-alppha)* mcs*kns^(alppha)*(1-tauns)*(1+tauds)^(el)) ...
    + (ydns*(1-omgs)*gsss*agds)/(1+tauds)^(eg) + delta*kns ...
    + ((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun))/((X(2)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)) ...
    * (zeta/zetas) * ( (X(1)*mun*oms*acm)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taum)^(el))...
    + (ydn*omgs*gss*agm)/(1+taum)^(eg) )) ,...
    X(3) - ((X(1) - muc*nuc)^(-ec)),...
    X(3)*(1+pi)*nubeta - (X(4) * (1+pis)* nubetas)];

% type 4 - modified
% eq = @(X) [ ( (X(1) *mun*(1-om)*acd)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taud)^(el)) ...
%     + (ydn*(1-omg)*gss*agd)/((1+taud)^(eg)) + delta * kn ...
%     + ((X(1)- muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns))/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)) ...
%     * (zetas/zeta) * ( (X(1)*mun*om*acms)/((X(1)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)*(1+taums)^(el)) ...
%     + (gsss*ydns*omg*agms)/(1+taums)^(eg) )) ...
%       - ( (X(1)*mun*(1-oms)*acds)/((X(1)-muc*nucs)^(-ec)*(1-alppha)* mcs*kns^(alppha)*(1-tauns)*(1+tauds)^(el)) ...
%     + (ydns*(1-omgs)*gsss*agds)/(1+tauds)^(eg) + delta*kns ...
%     + ((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun))/((X(1)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)) ...
%     * (zeta/zetas) * ( (X(1)*mun*oms*acm)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taum)^(el))...
%     + (ydn*omgs*gss*agm)/(1+taum)^(eg) )) ]

% type 2 
% eq = @(X) [ ydn - ( (X(1) *mun*(1-om)*acd)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taud)^(el)) ...
%     + (ydn*(1-omg)*gss*agd)/((1+taud)^(eg)) + delta * kn ...
%     + ((X(2)- muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns))/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)) ...
%     * (zetas/zeta) * ( (X(2)*mun*om*acms)/((X(2)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)*(1+taums)^(el)) ...
%     + (gsss*ydns*omg*agms)/(1+taums)^(eg) )) ,...
%     ydns - ( (X(2)*mun*(1-oms)*acds)/((X(2)-muc*nucs)^(-ec)*(1-alppha)* mcs*kns^(alppha)*(1-tauns)*(1+tauds)^(el)) ...
%     + (ydns*(1-omgs)*gsss*agds)/(1+tauds)^(eg) + delta*kns ...
%     + ((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun))/((X(2)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)) ...
%     * (zeta/zetas) * ( (X(1)*mun*oms*acm)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taum)^(el))...
%     + (ydn*omgs*gss*agm)/(1+taum)^(eg) )) ,...
%     X(3) - ((X(1) - muc*nuc)^(-ec)),...
%     X(4) - ((X(2) - muc*nuc)^(-ec))];

% type 3
% eq = @(X) [ ydn - ( (X(1) *mun*(1-om)*acd)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taud)^(el)) ...
%     + (ydn*(1-omg)*gss*agd)/((1+taud)^(eg)) + delta * kn ...
%     + ((X(2)- muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns))/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)) ...
%     * (zetas/zeta) * ( (X(2)*mun*om*acms)/((X(2)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)*(1+taums)^(el)) ...
%     + (gsss*ydns*omg*agms)/(1+taums)^(eg) )) ,...
%     ydns - ( (X(2)*mun*(1-oms)*acds)/((X(2)-muc*nucs)^(-ec)*(1-alppha)* mcs*kns^(alppha)*(1-tauns)*(1+tauds)^(el)) ...
%     + (ydns*(1-omgs)*gsss*agds)/(1+tauds)^(eg) + delta*kns ...
%     + ((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun))/((X(2)-muc*nucs)^(-ec)*(1-alppha)*mcs*kns^(alppha)*(1-tauns)) ...
%     * (zeta/zetas) * ( (X(1)*mun*oms*acm)/((X(1)-muc*nuc)^(-ec)*(1-alppha)*mc*kn^(alppha)*(1-taun)*(1+taum)^(el))...
%     + (ydn*omgs*gss*agm)/(1+taum)^(eg) )) ,...
%     X(1) - X(2),...
%     (X(3)*(1+pi)*nubeta)/(X(4) * (1+pis)* nubetas) - ( ((X(1) - muc*nuc)^(-ec) *(1+pi)*nubeta )/((X(2) - muc*nucs)^(-ec) *(1+pis)*nubetas) )];

tempX = fsolve(eq, [1.1, 1.1, 0.5, 0.5], options);
% tempX = fsolve(eq, [1.1], options);

x = tempX(1);                % total consumption in steady state (domestic)
xs = x;               % total consumption in steady state (foreign)
% lambda1 = tempX(3);    % marginal utility of consumption 'Blanchard calls it so' page 76.
% lambda1s = tempX(4);

lambda1 = (x - muc*nuc)^(-ec);
lambda1s = (xs - muc*nucs)^(-ec);

lambda2 = lambda1/nui;
lambda2s = lambda1s/nuis;

%  total habit adjusted consumption in steady state (domestic) 
c = x/(1-h);
%  total habit adjusted consumption in steady state (foreign) 
cs = xs/(1-h);

% variables which can be obtained by x directly
cd = ((x*(1-om))/(1+taud)^(el))* acd;                   % domestic consumption in steady state (domestic)
cm = ((x*om)/(1+taum)^(el))* acm;                       % imported consumption in steady state (domestic)
cds = (xs*(1-oms)*acds)/(1+tauds)^(el);              % domestic consumption in steady state (foreign)
cms = (xs*oms*acms)/(1+taums)^(el);                  % imported consumption in steady state (foreign)

% labor hours supply in steady state (domestic) ~ cant use it this way, it
% is dependent upon kn,
 n = ((x - muc * nuc)^(-ec)* (1-alppha) * mc* kn^(alppha)*(1-taun))/(mun* nun);
% labor supply in steady state (foreign)
 ns = ((xs-muc*nucs)^(-ec) * (1-alppha)* mcs * kns^(alppha)*(1-tauns))/(mun*nuns);
 
% capital can be obtained by the following:
k = n* kn;
ks = ns* kns;

% investment follows from capital
i = delta *k;
is = delta * ks;

% syms yd yds
% ydn = (cd/n) + gdn + (i/n) +(zetas/zeta)*(ns/n)*((cms/ns)+ gmns);
% ydns = (cds/ns) + gdns + (is/ns) + (zeta/zetas)*(n/ns)*((cm/n)+ gmn);
% [A,B] = equationsToMatrix([eq1, eq2], [yd, yds]);
% XF = linsolve(A,B);
% yd = double(XF(1,1));
% yds = double(XF(1,2));

yd      = ydn*n   ;             % domestic output supplied/demanded in steady state (domestic)
yds     = ydns*ns ;             % domestic output supplied/demanded in steady state (foreign)            


gd      = (gss*yd*(1-omg)*agd )/(1+taud)^(eg)  ;     % consumption on domestic goods by government at ss (domestic)
gm      = (gss*yd*omg*agm)/(1+taum)^(eg) ;             % consumption on imported goods by government at ss (domestic)
gds     = (gsss*yds*(1-omgs)*agds)/(1+tauds)^(eg) ;             % consumption on domestic goods by government at ss (foreign)
gms     = (gsss*yds*omgs*agms)/(1+taums)^(eg) ;             % consumption on imported goods by government at ss (foreign)
           


 

           

% numerator of reset price inflation in steady state (domestic)
f  = (nubeta * lambda1 *(1-tauyd)^(-ej)*mc*yd)/(1- xi*betta *((1+pid)/(1+pidss))^(ej));

% denominator of reset price inflation in steady state (domestic)
gpr = (nubeta * lambda1  *(1-tauyd)^(-ej)*yd)/(1- xi* betta * ((1+pid)/(1+pidss))^(ej-1));






g  = gss* yd      ;           % total governemnt spending in steady state (domestic)
gs = gsss * yds   ;           % total governemnt spending in steady state (foreign)



ydnat   = yd    ;             % natural rate of domestic output in steady state (domestic)
ydnats  = yds   ;             % natural rate of domestic output in steady state (foreign)










% numerator of reset price inflation in steady state (foreign)
fs = (nubetas* lambda1s*(1-tauydsss)^(-ej)*mcs*yds)/(1- xi*betta*((1+pids)/(1+pidss))^(ej));

% denominator of reset price inflation in steady state (foreign)
gprs = (nubetas * lambda1s * (1-tauydsss)^(-ej)*yds)/(1- xi*betta*((1+pids)/(1+pidss))^(ej-1));

% total output at steady state (cu)
ycu = zeta *yd + zetas * yds ;                  

% natural rate of total output at steady state (cu)
ycun = ycu ;                             

% output gap in steady state (cu)
ygapcu       = ycu - ycun ;            

% inflation in steady state (cu)
picu         = zeta * pi + zetas * pis ;   
picunat = pidss;
picugap = picu - picunat;


%%
%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state values of the endogenous variables of
%% the model.
%%

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%%
%% END OF THE SECOND MODEL INDEPENDENT BLOCK.