%compute the steady state of RBC Model 
function[ys,check]=erpt_steadystate(ys,exo)
global M_
 
%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for ii = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(ii,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(ii) ');']);         %    Get the value of parameter i.
end                                                           % End of the loop.  
check = 0;
%%


%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.
%%



 %Uses Analytical Solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta                =       0.9950;             %(2)   TIME-INVARIANT DISCOUNT FACTOR
L_ss                =       0.3000;             %(3)   STEADY STATE OF LABOUR HOURS
varepsilon          =       7.5000;             %(4)   ELASTICITY OF SUBSTITUTION BETWEEN THE DOMESTIC VARIETIES
epsilon             =       7.5000;             %(5)   ELASTICITY OF SUBSTITUTION BETWEEN THE LABOUR VARIETIES
delta               =       0.0250;             %(6)   CAPITAL DEPRECIATION RATE
nu_pi               =       1.5000;             %(7)   POLICY RATE RESPONSIVENESS TO INFLATION
nu_y                =       0.5000;             %(8)   POLICY RATE RESPONSIVENESS TO OUTPUT DEVIATION FROM THE STEADY STATE  
lambda              =       1.0000;             %(9)   RELATIVE SIZE OF DOMESTIC POPULATION
varsigma_pi         =       0.0100;             %(10)  DOMESTIC INFLATION TREND
gamma               =       0.6100;             %(11)  AVERAGE REAL EXCHANGE RATE
varsigma_y          =       1.5200;             %(12)  AVERAGE FOREIGN PER CAPITA INCOME RELATIVE TO THE DOMESTIC
%
omega               =       1.5000;             %(13)  CURVATURE OF THE UTILITY FUNCTION
varphi              =       3.0000;             %(14)  ELASTICITY OF LABOUR
chi                 =       0.5000;             %(15)  STRENGTH OF CONSUMPTION HABITS
eta                 =       1.5000;             %(16)  ELASTICITY OF SUBSTITUTION BETWEEN DOMESTIC AND FOREIGN GOODS
phi_p               =     120.0000;             %(17)  SIZE OF MENU COSTS
phi_w               =      50.0000;             %(18)  SIZE OF WAGE BARGAINING COSTS
phi_k               =       2.0000;             %(19)  SIZE OF INVESTMENT ADJUSTMENT COSTS
zeta_p              =   -1000.0000;             %(20)  DEGREE OF MENU COST ASYMMETRY
zeta_w              =    -350.0000;             %(21)  DEGREE OF WAGE BARGAINIGN COST ASYMMETRY
kappa               =       0.2300;             %(22)  CAPITAL SHARE OF INCOME NET OF COMMODITIES
xi                  =       0.2200;             %(23)  OPENNESS TO TRADE IN COMMODITIES
theta               =       0.4000;             %(24)  SHARE OF DISTRIBUTION SERVICES AT HOME
theta_star          =       0.4000;             %(25)  SHARE OF DISTRIBUTION SERVICES ABROAD
alpha               =       0.3000;             %(26)  DOMESTIC OPENNESS TO TRADE
alpha_star          =       0.1500;             %(27)  FOREIGN OPENNESS TO TRADE
%
rho_r               =       0.8100;             %(28)  PERSISTENCE OF NOMINAL INTEREST RATE
rho_a               =       0.7680;             %(29)  PERSISTENCE OF LABOUR PRODUCTIVITY
rho_y               =       0.8321;             %(30)  PERSISTENCE OF FOREIGN CONSUMPTION
rho_star            =       0.8000;             %(31)  PERSISTENCE OF FOREIGN PREFERENCES
sigma_lambda        =       0.0032;             %(32)  STANDARD ERROR OF MONETARY POLICY SHOCKS
sigma_a             =       0.0057;             %(33)  STANDARD ERROR OF PRODUCTIVITY SHOCKS
sigma_y             =       0.0064;             %(34)  STANDARD ERROR OF FOREIGN CONSUMPTION SHOCKS
sigma_star          =       0.0135;             %(35)  STANDARD ERROR OF FOREIGN PREFERENCE SHOCKS       
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PI=1+varsigma_pi;                                                                                       
PI_h=PI;                                                                                                
PI_tilde=1;                                                                                             
PI_W=1;                                                                                                 
q=gamma;                                                                                                
q_dot=1;                                                                                                
Q_dot=1;                                                                                                
r=1/beta;                                                                                               
R=r*PI;                                                                                                
%
DELTA_P=0;                                                                                              
dDELTA_P=0;                                                                                             
ddDELTA_P=phi_p;                                                                                        
DELTA_K=0;                                                                                              
dDELTA_K=0;                                                                                             
DELTA_W=0;                                                                                              
dDELTA_W=0;                                                                                             
TQ=1;                                                                                                   
x=(1-beta*(1-delta))/beta;                                                                              
LAMBDA=1;                                                                                               
A=1;                                                                                                   
%
x0=[1 1];
xis=fsolve(@nl_tot,x0);
XI_h=xis(1,1);
XI_f_star=xis(1,2);
%
THETA=XI_h;                                                                                             
THETA_star=XI_f_star;                                                                                   
s=(1+theta*THETA)/(1-theta);                                                                            
s_star=(1+theta_star*THETA_star)/(1-theta_star);                                                        
e=(gamma*XI_f_star)/(XI_h);                                                                             
MC=XI_h*((varepsilon-1)/varepsilon);                                                                    
%
kappa_M=xi;                                                                                             
kappa_K=(1-xi)*kappa;                                                                                   
kappa_N=1-kappa_M-kappa_K;                                                                              
%
GAMMA_k=(MC/lambda)*(kappa_K/x);                                                                        
GAMMA_m=(MC/lambda)*(kappa_M/q);                                                                        
GAMMA_y=((((GAMMA_k^kappa_K)*(GAMMA_m^kappa_M))*(lambda))^(1/kappa_N));                                 
GAMMA_w=(MC/lambda)*kappa_N*GAMMA_y;                                                                    
GAMMA_nx=(1-theta_star)*alpha_star*((XI_h/q)^(-eta))*varsigma_y-(1-theta)*alpha*(q*XI_f_star)^(-eta);                                                          
GAMMA_theta=((1-theta)*(1-alpha))*((s/XI_h)^eta);                                                       
GAMMA_c=(1-GAMMA_nx-delta*GAMMA_k*GAMMA_theta+gamma*GAMMA_m*GAMMA_theta)*(GAMMA_y/GAMMA_theta);         
%
w=GAMMA_w;                                                                                              
L=L_ss;                                                                                                 
N=lambda*L;                                                                                             
Y=(GAMMA_y*N)/GAMMA_theta;                                                                              
D=theta*Y;                                                                                              
Y_u=(1-theta)*Y;                                                                                        
Y_h=GAMMA_theta*Y;                                                                                      
K=GAMMA_k*Y_h;                                                                                          
M=GAMMA_m*Y_h;                                                                                          
I=delta*K;                                                                                              
C=GAMMA_c*N;                                                                                            
NX=GAMMA_nx*Y-q*M;                                                                                      
Y_star=varsigma_y*Y;                                                                                    
%
psi_ss=((GAMMA_w*(1-chi*beta)*(epsilon-1))/(varphi*epsilon))*L^(1-varphi);                              
%
PHI_V=(C*(1-chi)-psi_ss*(L^varphi))^(-omega);                                                           
UC=PHI_V*(1-chi*beta);                                                                                  
UC_star=UC;                                                                                             
UL=-psi_ss*varphi*(L^(varphi-1))*PHI_V;                                                                 
PHI_W_2=1;                                                                                             
PHI_W_1=-(epsilon/(epsilon-1))*(UL/UC);                                                                 
PHI_P_2=1;                                                                                              
PHI_P_1=(varepsilon/(varepsilon-1))*MC;                                                                 
F=Y-w*L-x*K;                                                                                            
OMEGA=-(NX)/(1-beta);                                                                                   
%
PHI_P_3=phi_p*(1+(beta/(PI^2)));                                                                        
%
PHI_Delta=PHI_P_3/(varepsilon-1);                                                                       
%
erpt=xi/(1+PHI_Delta);                                                                                  
erpt_star=1-erpt;                                                                                       
erpt_cpi_star=(1-theta_star)*alpha_star*((XI_h/q)^(1-eta))*erpt_star;                                   

%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state vZNues of the endogenous variables of
%% the model.
%%
NumberOfEndogenousVariables = M_.orig_endo_nbr;               % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for ii = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(ii,:));                      %    Get the name of endogenous variable i.                     
  eval(['ys(' int2str(ii) ') = ' varname ';']);                %    Get the steady state vZNue of this variable.
end                                                           % End of the loop.
%%