%%% steady state computation with numerically obtained parameters


function [ys,params,check] =ss_function(ys,exo,M_,options_)

	%global M_ 

	if isfield(M_,'param_nbr') == 1
        NumberOfParameters = M_.param_nbr;

        params=struct();

        for iter = 1:NumberOfParameters
            paramname = deblank(M_.param_names{iter,:});
            eval([ 'params.' paramname ' = M_.params(' int2str(iter) ');']);
        end
        check = 0;
	end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% find ss parameters by finding fixed point %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    guess_variables;
         
    guess_vars=Guess(:,1)';
    initial_guess=cell2mat(Guess(:,2))';

    N_guess=length(guess_vars);
    for j=1:N_guess
        guess_params.(guess_vars{j}) = initial_guess(j);
    end
    diff1=1;
    guess_old=initial_guess;
    iter_search=0;
    while diff1>10^-12&&iter_search<450
        iter_search=iter_search+1;
        for j=1:N_guess
            guess_old(j) = guess_params.(guess_vars{j}); % guess_old= [guess_params.kG_ss guess_params.yHstarVA_ss ...] 
        end 

        [upd] = update_ss(params,guess_params);          % upd = [ kG_ss_upd  yHstarVA_ss_upd ...]

        for j=1:N_guess
            guess_params.(guess_vars{j}) = upd(j);       % guess_params.kG_ss= kG_ss_upd ....
        end
        diff1  =((upd - guess_old)*(upd - guess_old)')^0.5;   % get root of the squared difference between updated and previous value of parameters       
    end
    
    DISP=0;

    
    if DISP==1
    disp(['stopped at ' num2str(iter_search) ' iterations, sum of squared differences: ' num2str(diff1^2)])
    end
    [ys,params_new]=get_steady(params,guess_params,M_);
    params=params_new;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% AUXILIARY FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [upd] = update_ss(params,guess_params)
%%% this function update the guessed parameters values for the steady state

	% Get the field names of the structure.
	fields.params = fieldnames(params, '-full');
	fields.guess_params = fieldnames(guess_params, '-full');
	% Find out how many there are - for our loop.
	numberOfFields.params = length(fields.params);
	numberOfFields.guess_params = length(fields.guess_params);

	for f = 1 : numberOfFields.params
		thisField = fields.params{f};
		commandLine = sprintf('%s = params.%s;', thisField, thisField);
		eval(commandLine);
	end

 	for f = 1 : numberOfFields.guess_params
 		thisField = fields.guess_params{f};
 		commandLine = sprintf('%s = guess_params.%s;', thisField, thisField);
 		eval(commandLine);
 	end

    %%%%%%%%%%%%% equations file %%%%%%%%%%%%
    steady_state_equations
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for f = 1 : numberOfFields.guess_params
 		thisField = fields.guess_params{f};
 		commandLine = sprintf('upd(%s) = %s_upd;', num2str(f), thisField);
 		eval(commandLine);
    end
        
    % Release temporary variables.
	clear('f', 'thisField', 'numberOfFields');
	clear('fields', 'commandLine');
end

function [ys,params_new] = get_steady(params,guess_params,M_)
%%% this function get the steady state values given a set of parameters

	%global M_  

	% Get the field names of the structure.
	fields.params = fieldnames(params, '-full');
	fields.guess_params = fieldnames(guess_params, '-full');
	% Find out how many there are.
	numberOfFields.params = length(fields.params);
	numberOfFields.guess_params = length(fields.guess_params);

	for f = 1 : numberOfFields.params
		thisField = fields.params{f};
		commandLine = sprintf('%s = params.%s;', thisField, thisField);
		eval(commandLine);
	end

 	for f = 1 : numberOfFields.guess_params
 		thisField = fields.guess_params{f};
 		commandLine = sprintf('%s = guess_params.%s;', thisField, thisField);
 		eval(commandLine);
 	end

	% Release temporary variables.
	clear('f', 'thisField', 'numberOfFields');
	clear('fields', 'commandLine');

    %%%%%%%%%%%%% equations file %%%%%%%%%%%%
    steady_state_equations
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for iter = 1:length(M_.params)
        eval([ 'M_.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);
        for iter = 1:NumberOfEndogenousVariables
            varname = deblank(M_.endo_names{iter,:});
            eval(['ys(' int2str(iter) ') = ' varname '_ss;']);
        end
    else
        check = 0;
    end
   
    x=zeros(M_.param_nbr,1);   %%% set exogenous variables to zero
    n_aux=M_.endo_nbr-M_.orig_endo_nbr;  %%% get number of aux variables
    ys=[ys;zeros(n_aux,1)];

    DISP=0;
    if DISP==1
    %%% compute residuals of the static model using the proposed ss
        %tic
        [RES,~,~] =eval([M_.fname '_static(ys, x, M_.params)']);
            for j=1:n_aux
                jj = M_.aux_vars(j).endo_index;
                ys(jj) = ys(jj)-RES(jj);
                [RES,~,~] =eval([M_.fname '_static(ys, x, M_.params)']);
            end    
        RES2=sum(RES.^2); %sum the residuals squares      
        disp(['Sum of the static model squared errors : ' num2str(RES2)])
        % toc 
    end
  params_new=M_.params;
end