% computes the steady state, including calibrating parameters
% calls PH_steadystate_sublevel which calculates steady state /given/ parameters

% if this is too slow, disable the subsidy calibration and/or Taylor approximations

function [ys,params,check] = PH_steadystate(ys, exe, M_, options_)

% we need oo_, which Dynare doesn't normally give to the steadystate file
% global oo_
% loc_oo = oo_; 

% Here we load the values of the deep parameters in a loop.

for ii = 1:M_.param_nbr                                 
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);         %    Get the value of parameter i.
end                                                        % End of the loop.  
check = 0;

% optimise over subsidy s while fixing n=0.5 with psys
% tsolve = @(x) cal_n_core(x, ys, M_);
% options = optimoptions('fsolve', 'Display', 'off');
% psy = fsolve(tsolve, psy, options);
% M_.params(find_string_index(M_.param_names, {'psy'})) = psy;

% given psy, optimise for subsidy s
tsolve = @(x) cal_s_core(x, ys, M_);
options = optimset('TolX', 1e-6, 'TolF', 1e-6);
s = fminsearch(tsolve, s, options); % set s to be utility-maximising 
M_.params(find_string_index(M_.param_names, {'s'})) = s;

%% now obtain the normalisation factor norc
[ys,~,~] = PH_steadystate_sublevel(ys, 1, M_);

% pull key variables 
for i = 1:length(M_.endo_names)
    eval([M_.endo_names{i} ' =  ys(' num2str(i) ');'])
end


% derive norc implied by U
norc = U / ((1-bet)*UU^gv+bet*U^gv)^(1/gv); 


WWbar = VM + VI + bbar*y;
vthe = (bbar*y*vthe_a)/WWbar;

% use SDF=bet and del1 = SDF*lamk/mu when u=1
del1 = mu / (bet*lamK);
%del1 = 26;

%% save everything done so far, before Taylor
% save calibrated parameters
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

for ii = 1:M_.orig_endo_nbr %auxiliary variables are set automatically
  eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
end

% save parameters to M_ for when we derive the Taylor approximations
M_.params = params;

%% Taylor approximations

% change options we need to force
change_options_ss;

% set decision rules 
dr = set_state_space(struct(),M_,options_);
dr.ys = ys;

UT1 = -66.166;
DT1 = -0.0553;

% loop Taylor approximations until converged
% for jj_temp = 1:2
%     
%     % generate decision rules using k_order_pert
%     orig_order = options_.order; % not sure why this is needed, but stochastic_solvers does it
%     [dr,info] = k_order_pert(dr,M_,options_);
%     options_.order = orig_order;
% 
%     % if error, print the error using Dynare's inbuilt error tool "print info"
%     if info(1)
%         disp(['Error triggered by steadystate file on loop' num2str(jj_temp)])
%         print_info(info, 0, options_)
%     end
%     
%     % run the approximation, outsourced to a subfile
%     obtain_Taylor_approx_second_order;
%     
%     % run steady_state subfile to update ups and delpie
%     [ys,~,check] = PH_steadystate_sublevel(ys, 1, M_);
%     dr.ys = ys;
% 
% end
% 



% run steady_state subfile to update ups and delpie
[ys,~, check] = PH_steadystate_sublevel(ys, 1, M_);