function [ys,params,check] = code_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = code_steadystate(ys,exo,M_,options_)
% computes the steady state for the code.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 <https://www.gnu.org/licenses/>.

% read out parameters to access them with their name

NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = 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

% INITIAL VALUES OF STEADY-STATE PARAMTERS
% betta=0.9434;        %domestic discount factor (6% rate of return).
% betta_star=0.9901;   %foreign discount rate (1% rate of return).
% delta=0.025;         %depreciation rate of private capital.
% delta_g=0.02;        %depreciation rate of public capital.
% alppha =0.333;       %non-oil out put elasticity to private capital
% alpa_o=0.7;          %Oil output elasticity to private capital.
% tau_o=0.27;          %Oil royalty rate.
% alpa_no=0.3;         %non-oil output elasticity to private capital
% lota_div=0.05;       %Dividend share of oil profit accrued to government.
% rho_r=0.95;          %interest rate smoothing in the Taylor rule   
% rho_rff=0.775;       %persistence in RFF accumlation process.
% zeta_GI=0.68;        %response of public investment to output.
% zeta_GC=0.30;        %response of public consumption to output.
% gammaa_GI =0.30;     %response of the public Investment to fiscal debt
% gammaa_GC =0.30;     %response of the public Consumption to fiscal debt
% gammaa_IOF=0.10;     %response of public invesment to oil revenue.
% gammaa_COF=0.20;     %response of public consumption to oil revenue
% phi=1.45;            %wage elasticity to hours worked
% phi_b=0.40;          %elasticity of fiscal bonds (debt)
% phi_I=0.20;          %elasticity of government investment
% phi_C=1.00;          %elasticity of government consumption
% phi_OR=-0.3;         %elasticity of Oil Revenue
% OMEGA=0.60;          %the upper bound of leverage ratio 0.54
% epi=9.00;            %elasticity of subsitution beteween different intermediate goods
% muo=0.5;             %the fraction of rule-of-thumb consumers in the economy
% PSI=0.16;            %non-oil output elasticity to public capital.
% sigmae=2.00;         %the inverse of the intertemporal elasticity of subsitution
% yobs=10;             %Observed (domestic) output
% cobs=10;             %Observed (domestic consumption
%set the parameter gammma1

% Steady State Values and Equations
    Z=1;                          %Technology (productivity) shock normalized to 1
    PX_d=1;                       %Domestic  inflation or CPI inflation
    PXstar=1;                     %foreign inflation
    PX_o=1;                       %Oil-price inflation
    PX_h=1;                       %Domestic inflation
    p_h=1;                        %domestic prices
    p_g=1;                        %Government purchase
    Pstar_o=1;                    %The world oil prices
    Ystar=1;   %500000;                 %We assume the steady state foreign output is zero 
    RER_SS=1;                     %If the law of one price holds, so the real exchange rate and relative prices equals to 1 in SS. 
    %KG_SS=3;                      %Capital (exogenous)*****Not sure the initial Public capital******
    Q_SS=1;                       %Steady-state Tobin Q (shadow price of capital)
    
    R_SS=1/betta;                 %Note that RSS=Rstar_SS=1/betta
    Rstar_SS=1/betta_star;        %Rstar_SS=R_SS;
    fxr_SS=betta_star-betta;
    Ko_SS=((1/betta-(1-delta))/(alpa_o*(1-tau_o)))^(1/(alpa_o-1));
    Yo_SS=(Ko_SS)^alpa_o;
    FDI_SS=delta*Ko_SS;
    RFF_SS=((tau_o+lota_div*(1-tau_o))*Yo_SS)/(1-rho_rff);
    OR_SS=(Rstar_SS-rho_rff)*RFF_SS*RER_SS;
    GI_SS=delta_g*KG_SS;
    b_SS=(((Yo_SS^(zeta_GI))*(OR_SS^(gammaa_IOF)))/(delta_g*KG_SS))^(1/gammaa_GI);
    GC_SS=((Yo_SS^(zeta_GC))*(OR_SS^(gammaa_COF))/(b_SS^gammaa_GC));
    T_SS=(((b_SS^(phi_b))*(GI_SS^(phi_I))*(GC_SS^(phi_C)))/(OR_SS^phi_OR)); %.XXX.causes colinear relationship
%     muo=1-((T_SS+OR_SS-GI_SS-GC_SS)/((R_SS-1)*b_SS)); %.XXX.causes colinear relationhip                      
    %T_SS =GI_SS+GC_SS+(1-muo)*(R_SS-1)*b_SS-OR_SS;
    Rkno_SS=(1/betta)-(1-delta)-(fxr_SS*OMEGA/betta);
    W_SS=(1-alppha)*(((KG_SS^PSI)*(alppha^alppha)*(epi-1))/(((Rkno_SS)^alppha)*epi))^(1/(1-alppha));
    N_SS=(W_SS)^(1/(phi-1));
    
% Y_SS=Yno_SS+(RER_SS*Yo_SS);  %the definition of Agggregate output (Y) as a sum of oil (Yo) and non-oil (Yno) output.     
% Kno_SS=(Rstar_SS*bstar_SS)/(Q_SS*PXstar); %New equation from steady state

%Kno_SS=(((RER_SS*Yo_SS))/((N_SS^(1-alppha))*(KG_SS^PSI)))^(1/alppha);  %Y_SS-
% Kno_SS=((Yno_SS)/((N_SS^(1-alppha))*(KG_SS^PSI)))^(1/alppha)
% Kno_SS=1;             %******NOT CLEAR NEED TO CALIBRATE**********
 KG_SS=((Y_SS-(RER_SS*Yo_SS))/((N_SS^(1-alppha))*(Kno_SS^alppha)))^(1/PSI); %Transformed non-oil capital equation.

   % KG_SS=(1/delta_g)*GI_SS;
    Yno_SS=((Kno_SS)^alppha)*(N_SS^(1-alppha))*(KG_SS^PSI); %since Z=1, this is the same as the notes
    %Kno_SS=((Y_SS-(RER_SS*Yo_SS))/((N_SS^(1-alppha))*(KG_SS^PSI)))^(1/alppha);
    Kno_SS=(1/delta)*I_SS;
    I_SS=delta*Kno_SS;
    bstar_SS=((OMEGA*Kno_SS))/Rstar_SS;
    NX_SS=(1-muo)*(Rstar_SS-1)*bstar_SS-(RER_SS*FDI_SS)-(Rstar_SS-rho_rff)*RFF_SS*RER_SS+(1-lota_div)*RER_SS*(1-tau_o)*Yo_SS;
    TS_SS=muo*((Rkno_SS-delta)*Kno_SS+b_SS*(R_SS-1)+bstar_SS*(1-Rstar_SS)+(1-((epi-1)/epi))*Yno_SS)+T_SS;
    TN_SS=(T_SS-(1-muo)*TS_SS)/muo;
    CN_SS=(W_SS*N_SS)-TN_SS;
    CS_SS=(W_SS*N_SS)-TN_SS+Rkno_SS+b_SS; %    CS_SS=CN_SS;
    C_SS=muo*CN_SS+(1-muo)*CS_SS;
    Y_SS=C_SS+(1-muo)*I_SS+GC_SS+GI_SS+NX_SS;
    
    
    
    R_SS=1/betta;
Rstar_SS=1; %%Rstar_SS=1/betta_star;
fxr_SS=betta_star-betta;
T_SS=(((b_SS^(phi_b))*(GI_SS^(phi_I))*(GC_SS^(phi_C)))/(OR_SS^phi_OR));
C_SS=muo*CN_SS+(1-muo)*CS_SS;
N_SS=(W_SS)^(1/(phi-1));



    
%  lambda=((1-betta*theta)*(1-theta))/theta; %definition of lambda   
%  lambda=(CS_SS-(N_SS^phi)/phi)^-sigmae;

%------------------------------------------------------------------------------------------------------
% SAVE STEADY-STATE VALUES OF VARIABLES
T= log(T_SS);       %T_SS = T;        
TS= log(TS_SS);     %TS_SS = TS;
Y= log(Y_SS);       %Y_SS=Y; 
Yo=log(Yo_SS);      %Yo_SS=Yo;  
Yno= log(Yno_SS);   %Yno_SS=Yno; 
Kno= log(Kno_SS);   %Kno_SS=Kno; 
Ko= log(Ko_SS);     %Ko_SS=Ko; 
b= log(b_SS);       %b_SS=b; 
bstar= log(bstar_SS); %bstar_SS=bstar; 
C= log(C_SS);       %C_SS=C; 
CS= log(CS_SS);     %CS_SS=CS; 
CN= log(CN_SS);     %CN_SS=CN; 
N= log(N_SS);       %N_SS=N;
% W= log(W_SS);       %W_SS=W;
NX= log(NX_SS);     %NX_SS=NX; 
R= log(R_SS);       %R_SS=R; 
Rstar= log(Rstar_SS);% Rstar_SS=Rstar; 
GI= log(GI_SS);     %GI_SS=GI; 
GC= log(GC_SS);     %GC_SS=GC; 
KG= log(KG_SS);     %KG_SS=KG; 
RER= log(RER_SS);   %RER_SS=RER; 
%OR= log(OR_SS);     %OR_SS=OR; 
I= log(I_SS);       %I_SS=I; 
fxr_star= log(fxr_SS);%fxr_SS=fxr;
RFF= log(RFF_SS);   %RFF_SS=RFF; 
FDI= log(FDI_SS);   %FDI_SS=FDI;
W= log(W_SS);       %W_SS=W; 
%%
%% 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.




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

if isfield(M_,'param_nbr') == 1
if isfield(M_,'orig_endo_nbr') == 1
  
NumberOfEndogenousVariables = M_.orig_endo_nbr;
else
NumberOfEndogenousVariables = M_.endo_nbr;
end
ys = zeros(NumberOfEndogenousVariables,1);
if options_.linear
    else
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
eval(['ys(' int2str(ii) ') = ' varname ';']);
end
end
else
ys=zeros(length(lgy_),1);
for ii = 1:length(lgy_)
    ys(ii) = eval(lgy_(ii,:));
end
check = 0;
end
