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 

a = omega*a1 + (1-omega)*a2 ;
b = omega*b1 + (1-omega)*b2 ;
%L_r	= 1 ;
%L_m 	= (1-omega)*L_r ;
%L_p	= 1 ;
i_s 	= (1/beta) - 1 ;
i_d 	= i_s ;
i_g 	= i_d ;
i_L	= (i_d - alpha_s*i_g)*(1 - alpha_c - alpha_s) ;

SC	= 1 ;
X_p 	= X_p_bar ;
inve	= 1 ;
macro	= 1 ;
shock_omss = 1 ;
shock_pds = 1 ;
pi_w	= 1 ; 		%% assuming pi_w to be in level which will log afterwards
pi	= 1 ;

P_v	= 1 ;
MC_v	= (psi1 - 1)/psi1 ;
A_v	= 1 ;
A 	= 1 ;
A_m 	= 1 ;


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


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

%[FGP,fval,exitflag] = fzero(@(FGP) (FGP + pc*((alpha1^(alpha1/(1-alpha1-beta1)))*(beta1^(beta1/(1-alpha1-beta1)))*(FGP^((alpha1+beta1)/(1-alpha1-beta1)))*((w_p)^((-beta1)/(1-alpha1-beta1)))*(r^((-alpha1)/(1-alpha1-beta1)))) - (2*beta)/((X_p_bar)*(1- beta + beta*del_p)) ), 1,options) ; 

MC_g	= (1/A)*((w_p/(1-alpha1))^(1-alpha1))*((r/alpha1)^alpha1) ;
FGP 	= MC_g ;
MSP 	= FGP ;
Y_g_D	= ( (2*beta)/((X_p_bar)*(1- beta + beta*del_p)) - FGP )/2*pc ; 
Y_g_FGP	= Y_g_D/(1 - m) ;

%FGP	= (2*beta)/((X_p_bar)*(1- beta + beta*del_p)) - 2*pc*Y_g_D ;

%Y_g_MSP	= Y_g_FGP ;


Y_g_S	= Y_g_D - del_p*X_p_bar ;

%L_g_MSP	= (1-alpha1)*(MC_g)*Y_g_MSP/w_p ;
%K_g_MSP	= alpha1*(MC_g)*Y_g_MSP/r ;

K_g_FGP	= alpha1*(MC_g)*Y_g_FGP/r ;
L_g_FGP	= (1-alpha1)*(MC_g)*Y_g_FGP/w_p ;

Y_g_MSP	= (proc/(1-proc))*((1-m)/m)*Y_g_FGP ;

K_g_MSP	= alpha1*(MC_g)*Y_g_MSP/r ;
L_g_MSP	= (1-alpha1)*(MC_g)*Y_g_MSP/w_p ;

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

L_g	= m*L_g_MSP + (1-m)*L_g_FGP ;
K_g	= m*K_g_MSP + (1-m)*K_g_MSP ;
Y_g	= m*Y_g_MSP + (1-m)*Y_g_MSP ;

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

%% Now the Govt. sector

PDS 	= offtake*omega ;
X_g 	= (m*Y_g_MSP - PDS + X_g_bar)/(1+del_g) ;
OMSS	= (m*Y_g_MSP - PDS - del_g*X_g_bar)/(1+del_g) ; 

G 	= (1+mu)*MSP*m*Y_g_MSP ;
tax	= G/(1 - lambda*i_g) ;
B_s	= lambda*tax ;
Pss_Css	= (1+mu)*MSP*m*Y_g_MSP/(tau*(1-a-b)*(1 + lambda*i_g)) ;

P_g	= a*Pss_Css/(Y_g_S + OMSS + omega*PDS) ;
Y_v	= b*Pss_Css ;

L_v	= (1-alpha2)*MC_v*Y_v/w_p ;
K_v	= alpha2*MC_v*Y_v/r ;

L_p	= (L_v + L_g)/omega ;




%% Industrial sector 
L_r	= L_p ;
L_m	= (1-omega)*L_r ;
w_r	= (psi_w/(psi_w - 1))*(1/(1-beta*theta_w))*w_p ;
w_star_r = w_r ;
w_m 	= w_r ;
MC_m	= ((w_r/(1-alpha3))^(1-alpha3))*((r/alpha3)^alpha3) ;
P_m 	= (1/(1-tau))*(psi2/(psi2-1))*(1/(1-beta*theta))*MC_m ;
P_star_m = P_m ;
Y_m	= tax/(tau*P_m) ; 	%% since tax = tau*P_m*Y_m

%L_m	= (1-alpha3)*MC_m*Y_m/w_m ; %%Problem: there are 2 ways to dtermine L_m 
K_m	= (alpha3)*MC_m*Y_m/r ;



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

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

%% consumption and savings of HH

C	= Pss_Css/P ;
C_p	= C
C_r	= C ;

% Now let 

%S	= (P_r*C_r - w_r*L_r)/i_s ; %% only rich HH save 
%B_d	= (i_g/(i_d + i_g))*(1-omega)*S ;
%D 	= (i_d/(i_d + i_g))*(1-omega)*S ;


S 	= ((i_d + i_g)*B_s)*(1/(alpha_s*i_d + i_g))*(1/(1-omega)) ; 	%check where to add 1-omega
B_d	= (i_g/(i_d + i_g))*(1-omega)*S ;
D 	= (i_d/(i_d + i_g))*(1-omega)*S ;

C_g_p 	= a1*(P_p/P_g)*C_p - PDS ;
C_v_p 	= b1*(P_p/P_v)*C_p ; 	
C_m_p 	= (1-a1-b1)*(P_p/P_m)*C_p ;

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

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



%% Banking and capital producing firm

K 	= K_g + K_v + K_m ;
I	= del_k*K ;
loan	= I ;
Z	= (1/del_b)*(i_L*loan + i_g*alpha_s*D - i_d*D) ;
gpft	= del_b*Z ;





%% 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 ';']);

eval(['ys(' int2str(ii) ') = log(' varname ');']);
end
