function [residuals] = BBCDLmodelcomp_ramsey_c1pid_c2pid_endoss_main(instr1,instr2,M_,ys,keep_variables,exogenous_variables)

% except for function names, this file does not depend on the model 

global ys_conditional info_lmss_mat ys_guess
global iteration

% to reduce the complexity of the problem we rewrite the steady state
% conditions so that we only have to provide initial guesses for a subset
% of the endogenous variables


%% Definitions and assignments
for index = 1:length(keep_variables)
    keep_variables_pos(index) = strmatch(keep_variables{index},M_.endo_names,'exact');
end
M_.endo_names_reduced = M_.endo_names(keep_variables_pos,:);


%%% CHECK THIS PART MAY DEPENND ON MODEL%%%
keep_nonlm_pos =[];
for idxxx=1:1:size(M_.endo_names,1)
    if isempty(strmatch('lmult',M_.endo_names(idxxx,1:5),'exact'))
        keep_nonlm_pos = [keep_nonlm_pos;idxxx];
    end
end
M_.endo_names_conditional = M_.endo_names(keep_nonlm_pos,:);

% for open-loop Nash
% keep_nonlm_pos =[];
% for idxxx=1:1:size(M_.endo_names,1)
%     if isempty(strmatch('lmult',M_.endo_names(idxxx,3:7),'exact'))
%         keep_nonlm_pos = [keep_nonlm_pos;idxxx];
%     end
% end
% M_.endo_names_conditional = M_.endo_names(keep_nonlm_pos,:);



%% obtain solution for keep_variables
variablesguess = ys(keep_variables_pos,:);
algo = char('active-set', 'trust-region-reflective', 'interior-point', 'interior-point-convex', 'levenberg-marquardt', 'trust-region-dogleg', 'lm-line-search','sqp');
options = optimset('Display','Iter','MaxFunEvals',1e4,'MaxIter',1e4,'TolFun',1e-9,'TolX',1e-9,'Algorithm',algo(5,:));
[ys_reduced,fval,exitflag]=fsolve(@(variablesguess) BBCDLmodelcomp_ramsey_c1pid_c2pid_endoss_reduced1(variablesguess,instr1,instr2,exogenous_variables,M_),variablesguess,options);



if exitflag < 1
    residuals = [NaN;NaN];
    return
end
    

% assign new steady state values to all (structural) endogenous variables
for idxxx=1:1:size(M_.endo_names_conditional,1)
    define_string=[deblank(M_.endo_names_conditional(idxxx,:)),'_ss =',sprintf('%20.20g',ys_conditional(idxxx))];
    evalc(define_string);
end

if isnan(M_.params(:,:))==0
    %read parameter values from M_.params and assign to parameters names.
    %this overwrites parameters coming from params.m if estimation is done.
    num_para=size(M_.params,1);
    for idxxx=1:1:num_para
        define_string=[M_.param_names(idxxx,:),' =',sprintf('%20.20g',M_.params(idxxx))];
        evalc(define_string);
    end
end

%% initialize first order conditions with respect to the endogenous variables to compute lagrange multipliers
% the first order conditions are linear in the lagrange multipliers
% lmss_vec - lmss_mat *lmss_val = 0
% the system is overdetermined with lmss_mat of dimension N by (N-2) and 
% lmss_val of dimension (N-2) by 1.

% If the steady state inflation rate changed and we want to make indexation
% consistent with this number, we need to alter the parameter file
% accordingly:
c1pistar = c1pid_ss;
c2pistar = c2pid_ss;

BBCDLmodelcomp_ramsey_c1pid_c2pid_lmss


% check rank condition for coefficient matrix
lmss_mat_rank = rank(lmss_mat);
[n1_lmss,n2_lmss]= size(lmss_mat);

% find a submatrix of full rank
% first create a matrix of full row rank
% ind1_removed records the equations that are removed from the system at
% this stage (need to check them later)
[q,r,e1]=qr(lmss_mat');
[junk,ind1]=max(e1);
ind1_sorted = sort(ind1(1:lmss_mat_rank));
ind1_removed = sort(ind1(lmss_mat_rank+1:end)); 
lmss_mat_fullrowrank = lmss_mat(ind1_sorted,:);
lmss_mat_fullrowrank_prime = lmss_mat_fullrowrank';

% second create a matrix of full row and column rank
% if lmss_mat_rank < N-2, we need to reduce lmss_mat_fullrowrank_prime in 
% dimension to obtain a matrix of full row and column rank
% ind2_removed records the columns that are removed from the system
% if ind2_removed > 0, not all lagrange multipliers are determined
% (multiple solutions from lm)
[q,r,e2]=qr(lmss_mat_fullrowrank_prime');
[junk,ind2]=max(e2);
ind2_sorted = sort(ind2(1:lmss_mat_rank));
ind2_removed = sort(ind2(lmss_mat_rank+1:end)); 
lmss_mat_fullcolumnrank_prime = lmss_mat_fullrowrank_prime(ind2_sorted,:);
lmss_mat_fullrank = lmss_mat_fullcolumnrank_prime';
lmss_mat_fullrank_removedrow = lmss_mat_fullrowrank_prime(ind2_removed,:)';


lmss_vals_fullrank = zeros(n2_lmss,1);
if ~isempty(ind2_removed)

    if iteration == 1
        lmss_vals_fullrank(ind2_removed,1) = zeros(length(ind2_removed),1);
    elseif iteration == 2
        lmss_vals_fullrank(ind2_removed,1) = zeros(length(ind2_removed),1)+1.9;
    end
end
lmss_vals_fullrank(ind2_sorted,1) = -inv(lmss_mat_fullrank)*(lmss_vec(ind1_sorted,1)...
                                        +lmss_mat_fullrank_removedrow*lmss_vals_fullrank(ind2_removed,1));

% unused equations and undetermined lagrange multipliers
residuals = lmss_mat(ind1_removed,:)*lmss_vals_fullrank - lmss_vec(ind1_removed,1);

%% store current values in ys_conditional
for idxxx=1:1:size(M_.endo_names_conditional,1)
    define_string=[deblank(M_.endo_names_conditional(idxxx,:)),' =',sprintf('%20.20g',ys_conditional(idxxx))];
    evalc(define_string);
end

for idxxx=1:1:size(lmss_vals,1)
    define_string=['lmult',num2str(idxxx),' =',sprintf('%20.20g',lmss_vals_fullrank(idxxx))];
    evalc(define_string);
end    
    
for indxxx = 1:1:size(M_.endo_names,1)
    ys_guess(indxxx,1) = eval(M_.endo_names(indxxx,:));
end

% other information
lmss_mat_rank = rank(lmss_mat);
lmss_mat_fullrank_rank = rank(lmss_mat_fullrank);
[numeq,numlm] = size(lmss_mat);
numallzerorows = sum(sum(abs(lmss_mat),2)==0);
info_lmss_mat = [lmss_mat_rank,lmss_mat_fullrank_rank,numeq,numlm,numallzerorows];


