%Steady state computation file for Kedia(2021) base model
%NOTE: Not tested with Dynare versions after 4.5.7

function [ys,params,check] = K2021_BASE_EXP_steadystate(ys,exe,M_,options_)
%% 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.
  eval([ 'paramstemp.' paramname '= M_.params(' int2str(ii) ');']);
end                                                           % End of the loop.  
check = 0;
%%
%% END OF THE FIRST MODEL INDEPENDENT BLOCK.


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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
starting.K0 = K0;
starting.L0 = L0;
starting.chi0 =  chi0;            
starting.lambda0  =  lambda0;     
starting.omega0 = omega0;

%Setting target moments
moments.L_mom    =   L_mom;
moments.RkmR_mom =   RkmR_mom;
moments.phi_mom  =   phi_mom;
switch_print    =   'no';
switches.switch_print    =  switch_print;  %'test' for too many details, 'details', 'yes', 'no'

%Calibrating the parameters to hit the moments
switch switch_print
    case 'details'
        options     =   optimset('Display','iter');
    case 'test'
        options     =   optimset('Display','iter');
    otherwise
        options     =   optimset('Display','off');
end;

[XX,diff,exitf]     =   fsolve(@f_mom,[omega0 lambda0 chi0],options,paramstemp,starting,moments,switches);


omega  =   XX(1);
lambda =   XX(2);
chi    =   XX(3);


paramstemp.omega    =   omega;
paramstemp.lambda   =   lambda;
paramstemp.chi      =   chi;

%Calculating the equilibrium K and N values given the parameters
[XX_ss,fval,exitf]  =   fsolve(@f_KL,[K0 L0],options,paramstemp);

K    =   log(XX_ss(1));
L    =   log(XX_ss(2));

switch switches.switch_print
    case 'test'
        options     =   optimset('Display','iter');
    otherwise
        options     =   optimset('Display','off');
end;

%Miscellaneous values
Pm   =  log((epsilon-1)/epsilon);
X    =  log(1/((epsilon-1)/epsilon));
D    =  log(1);
Q    =  log(1);
delta= log(deltai);

%Macrovariables
Y   =  log(exp(K)^alfa*exp(L)^(1-alfa)); 
Y_ss= Y;
Ym  =  Y;
I   =  log(deltai*exp(K));
I_ss= exp(I);
In  =  0;
G   =  log(G_over_Y*exp(Y));
G_ss= exp(G);
C   =  log(exp(Y)-exp(I)-exp(G)); 
varrho  =   log((1-betta*hh)*((1-hh)*exp(C))^(-sig));  
Lambda  =  log(1);
R       =  log(1/betta);
Rk      =  log(exp(Pm)*alfa*exp(Y)/exp(K)+1-deltai);
Rs      =  Rk;
RkmR    =  log(exp(Rk)-exp(R));
w       =  log(exp(Pm)*(1-alfa)*exp(Y)/exp(L));  
VMPK    =  log(exp(Pm)*alfa*exp(Y)/exp(K));      
U       =  log(1);

%Variable capacity utilization parameters
b       =    exp(Pm)*alfa*exp(Y)/exp(K);  
delta_c=exp(delta)-b/(1+zetta);

%Financial Intermediaries
aa     =   lambda*betta*theta*exp(RkmR);
bb     =   -(1-theta)*(lambda-betta*exp(RkmR));
cc     =   (1-theta);
phi    =   log((-bb-sqrt(bb^2-4*aa*cc))/(2*aa));
z      =   log(exp(RkmR)*exp(phi)+exp(R));
x      =   z;         
nu     =   log(((1-theta)*betta*exp(RkmR))/(1-betta*theta*exp(x)));
eta    =   log((1-theta)/(1-betta*theta*exp(z)));
N      =   log(omega*exp(K)/(1-theta*(exp(RkmR)*exp(phi)+exp(R))));
Ne     =   log(theta*exp(z)*exp(N));
Nn     =   log(omega*exp(K)); 
prem   =   log(exp(Rk)/exp(R));

%Pricing variables
F   =   log(exp(Y)*exp(Pm)/(1-betta*gam));
Z   =   log(exp(Y)/(1-betta*gam));
infl  =   0;
inflstar =0;
i   =   R;
 
Y_OBS = ctrend;
I_OBS = ctrend;
C_OBS = ctrend;
PI_OBS = infl*100 + pi_star;
i_OBS = i*100 + i_star;
N_OBS = 0;
Rs_OBS = beta_0 + beta_1*2*Rs;

%Shocks
a=0;
ksi=0; 
g=0;
                             
e_a     =   0;
e_ksi   =   0;
e_g     =   0;
e_Ne    =   0;
e_i     =   0;
e_conf  =   0;
e_meas  =   0;

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

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 value of this variable.
end
 