function [ys,check] = MSP_steadystate(ys,exo)   

% function [ys,check] = MSP_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;




%% Enter model equations here


options=optimset(); % set options for numerical solver



% the steady state computation follows 

exp(L_r)	= 1 - omega ;
exp(L_m) 	= exp(L_r) ;
exp(L_p)	= omega ;
exp(i_s) 	= (1/beta) - 1 ;
exp(i_d) 	= exp(i_s) ;
exp(i_g) 	= exp(i_d) ;

exp(P_v)	= 1 ;
exp(MC_v)	= (psi1 - 1)/psi1 ;
exp(A_v)	= 1 ;
exp(A) 		= 1 ;
exp(A_m) 	= 1 ;
exp(i)		= ((1-alpha_s)*exp(i_s))/(1-alpha_c-alpha_s) ;
exp(i_L) 	= exp(i) ;

exp(r) 		= 1 + exp(i_L) ;
exp(w_v) 	= (1-alpha2)*(((psi1-1)/psi1)^(1/(1-alpha2)))*((alpha2/exp(r))^(alpha2/(1-alpha2))) ;



exp(m) 		= 0.5 ;

%exp(FGP)	= fzero(exp(FGP) + pc*(((alpha1)^(alpha1/(1-alpha1-beta1))*((beta1)^(beta1/(1-alpha1-beta1)))*((exp(FGP))^((alpha1+beta1)/(1-alpha1-beta1)))*((exp(w_g))^((-beta1)/(1-alpha1-beta1)))*((exp(r))^((-alpha1)/(1-alpha1-beta1)))) - (2*beta)/((X_p_bar)*(1- beta + beta*del_p)), 1) ;
%fun 		= exp(FGP) + pc*(((alpha1)^(alpha1/(1-alpha1-beta1))*((beta1)^(beta1/(1-alpha1-beta1)))*((exp(FGP))^((alpha1+beta1)/(1-alpha1-beta1)))*((exp(w_g))^((-beta1)/(1-alpha1-beta1)))*((exp(r))^((-alpha1)/(1-alpha1-beta1)))) - (2*beta)/((X_p_bar)*(1- beta + beta*del_p)) ;
%exp(FGP)	= fzero(fun, 1) ;

%exp(FGP)	= fzero(exp(FGP) + pc*(((alpha1)^(alpha1/(1-alpha1-beta1))*((beta1)^(beta1/(1-alpha1-beta1)))*((exp(FGP))^((alpha1+beta1)/(1-alpha1-beta1)))*((exp(w_g))^((-beta1)/(1-alpha1-beta1)))*((exp(r))^((-alpha1)/(1-alpha1-beta1)))) - (2*beta)/((X_p_bar)*(1- beta + beta*del_p)), 1) ;

%exp(Y_g_FGP)	= exp(m)*((alpha1)^(alpha1/(1-alpha1-beta1))*((beta1)^(beta1/(1-alpha1-beta1)))*((exp(FGP))^((alpha1+beta1)/(1-alpha1-beta1)))*((exp(w_g))^((-beta1)/(1-alpha1-beta1)))*((exp(r))^((-alpha1)/(1-alpha1-beta1))) ;

exp(Y_g_MSP) 	= exp(m)*((alpha1^(alpha1/(1-alpha1-beta1)))*(beta1^(beta1/(1-alpha1-beta1)))*(exp(FGP)^((alpha1+beta1)/(1-alpha1-beta1)))*(exp(w_p)^((-beta1)/(1-alpha1-beta1)))*(exp(r)^((-alpha1)/(1-alpha1-beta1)))) ;

exp(FGP)	= (2*beta)/((X_p_bar)*(1- beta + beta*del_p)) - 2*pc*exp(Y_g_MSP) ;

exp(Y_g_D)	= 0.5*exp(Y_g_FGP) ;
exp(Y_g_MSP)	= exp(Y_g_FGP) ;

%exp(Y_g_FGP)	= 0.5*exp(Y_g) ;

exp(Y_g_S)	= exp(Y_g_D) - del_p*X_p_bar ;

% exp(L_g)	= ((alpha1)^(alpha1/(1-alpha1-beta1))*((beta1)^((1-alpha1)/(1-alpha1-beta1)))*((exp(FGP))^((1)/(1-alpha1-beta1)))*((exp(w_g))^((-(1-alpha1))/(1-alpha1-beta1)))*((exp(r))^((-alpha1)/(1-alpha1-beta1))) ;

exp(L_g) = ((alpha1^alpha1)*(beta1^(1-alpha1))*(exp(FGP))*(exp(w_p)^(alpha1-1))*(exp(r)^(-alpha1)))^(1/(1-alpha1-beta1)) ;

% exp(K_g)	= ((alpha1)^((1-beta1)/(1-alpha1-beta1))*((beta1)^(beta1/(1-alpha1-beta1)))*((exp(FGP))^((1)/(1-alpha1-beta1)))*((exp(w_g))^((-beta1)/(1-alpha1-beta1)))*((exp(r))^((-(1-beta1))/(1-alpha1-beta1))) ;

exp(K_g) = ((alpha1^(1-beta1))*(beta1^(beta1))*(exp(FGP))*(exp(w_p)^(-beta1))*(exp(r)^(beta1-1)))^(1/(1-alpha1-beta1)) ;

exp(K_g_MSP)	= 0.5*exp(L_g) ;
exp(K_g_FGP)	= 0.5*exp(L_g) ;
exp(L_g_MSP) 	= 0.5*exp(L_g) ;
exp(L_g_FGP) 	= 0.5*exp(L_g) ;
exp(MSP) 	= exp(FGP) ;


exp(L_v) 	= exp(L_p) - exp(L_g) ;
exp(Y_v)	= (exp(w_v)*exp(L_v))/((1-alpha2)*exp(MC_v)) ; 
exp(K_v)	= alpha2*exp(MC_v)*(exp(Y_v)/exp(r)) ;

%%industrial sector 
exp(w_r)	= (exp(Y_v) - omega*b1*exp(w_p))/((1-omega)*b2*(1-beta*theta_w)*((psi_w-1)/psi_w)) ;
exp(w_m) 	= exp(w_r) ;
exp(MC_m)	= ((exp(w_r)/(1-alpha3))^(1-alpha3))*((exp(r)/alpha3)^alpha3) ;
exp(P_m) 	= (1/(1-tau))*(psi2/(psi2-1))*(1/(1-beta*theta))*exp(MC_m) ;
exp(Y_m)	= (exp(w_m)*exp(L_r))/((1-alpha3)*exp(MC_m)) ; 
exp(K_m)	= exp(Y_m)^(1/alpha3) ;

%% Now P_g remaining 

exp(OMSS)	= 0;
exp(P_g) 	= (omega*a1*exp(w_p) + (1-omega)*a2*((psi_w - 1)/psi_w)*(1-beta*theta_w)*exp(w_r))/exp(Y_g_S) ;

%% now we have all prices, so we compute all price indices

exp(P_p) 	= ((exp(P_g))^a1)*((exp(P_v))^b1)*((exp(P_m))^(1-a1-b1)) ;
exp(P_r) 	= ((exp(P_g))^a2)*((exp(P_v))^b2)*((exp(P_m))^(1-a2-b2)) ;
exp(P) 	= (exp(P_p)^omega)*(exp(P_r)^(1-omega)) ;

%% consumption and savings of HH

exp(C_p)	= exp(w_p)/exp(P_p) ;
exp(C_r)	= ((psi_w-1)/psi_w)*(1-beta*theta_w)*(exp(w_r)/exp(P_r)) ;

exp(S)		= (exp(P_r)*exp(C_r) - exp(w_r)*exp(L_r))/exp(i_s) ; %only rich HH save 
exp(B_d)	= (exp(i_g)/(exp(i_d) + exp(i_g)))*exp(S) ;  	%check for extra bracket
exp(D) 	= (exp(i_d)/(exp(i_d) + exp(i_g)))*exp(S) ;


exp(shock_pds)	= 1 ;
exp(PDS) 	= offtake*omega + exp(shock_pds) ;   %reconsider shock part
exp(C_g_p) 	= a1*(exp(P_p)/exp(P_g))*exp(C_p) - exp(PDS) ;
exp(C_v_p) 	= b1*(exp(P_p)/exp(P_v))*exp(C_p) ; 	
exp(C_m_p) 	= (1-a1-b1)*(exp(P_p)/exp(P_m))*exp(C_p) ;

exp(C_g_r) 	= a1*(exp(P_r)/exp(P_g))*exp(C_r) ;
exp(C_v_r) 	= b1*(exp(P_r)/exp(P_v))*exp(C_r) ; 	
exp(C_m_r) 	= (1-a1-b1)*(exp(P_r)/exp(P_m))*exp(C_r) ; 

exp(C_g) 	= omega*exp(C_g_p) + (1-omega)*exp(C_g_r) ;
exp(C_v) 	= omega*exp(C_v_p) + (1-omega)*exp(C_v_r) ;
exp(C_m) 	= omega*exp(C_m_p) + (1-omega)*exp(C_m_r) ;



%% Banking and capital producing firm

exp(K) 		= exp(K_g) + exp(K_v) + exp(K_m) ;
exp(I)		= del_k*exp(K) ;
exp(loan)	= exp(I) ;
exp(Z)		= kb*exp(loan) ;
exp(B_IB)	= (1-kb)*exp(loan) - (1 - alpha_c - alpha_s)*exp(D) ;
exp(gpft)	= del_b*exp(Z) ;

%% Now the Govt. sector

exp(tax) 	= tau*exp(P_m)*exp(C_m) ;
exp(X_g_bar)	= (exp(m)*exp(Y_g_MSP)  - exp(OMSS) - exp(PDS))/del_g ;			%check what is X_g_bar?
 

exp(B_s)	= (exp(tax) - exp(MSP)*exp(m)*exp(Y_g) - (proc*(exp(m)*exp(Y_g))^2) - stor*((exp(X_g_bar))^2) - dist*((exp(PDS))^2))/exp(i_g) ;




%% 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
