function [fval] = BBCDLmodelcomp_ramsey_c1pid_c2pid_endoss_reduced(variablesguess,instr1,instr2,exogenous_variables,M_)
 
% this file is customized to the model under consideration to speed up
% computations

global ys_conditional

%% assign variables and parameter values 
num_variablesguess=size(variablesguess,1);
for idxxx=1:1:num_variablesguess
    define_string=[M_.endo_names_reduced(idxxx,:),'=',sprintf('%20.20g',variablesguess(idxxx))];
    evalc(define_string);
end

%exogenous_variables = cellstr({'c1z';'c1taupvar';'c1shockgy';'c2z';'c2taupvar';'c2shockgy'});
num_exvar=size(exogenous_variables,1);
for idxxx=1:1:num_exvar
    define_string=[char(exogenous_variables(idxxx,:)),'=',sprintf('%20.20g',0)];
    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

%% compute other endogenous variables conditional on variablesguess

% needed to evaluate fval
c1pid   = instr1;
c2pid   = instr2;
% %make inflation indexation consistent with steady state inflation rate
% c1pistar = c1pid;
% 
% c2pistar = c2pid;
% 
c1rpcpd = log(1/(c1omegac*exp(c1c)/exp(c1cd))^(c1rhoc/(1+c1rhoc))); 
c2rpcpd = log(1/(c2omegac*(exp(c2c))/(exp(c2cd)))^(c2rhoc/(1+c2rhoc))); 
c1rer   = log((c1rerss*c1lambdaqss/((exp(c1c))^(-c1sigma))*c1rpcpdss)/(c2lambdaqss/((exp(c2c))^(-c2sigma))*c2rpcpdss)); 
c1mc    = log(c1omegamc*exp(c1c)/(((exp(c1rer))/(exp(c2rpcpd)))^((1+c1rhoc)/c1rhoc))); 
c2mc    = log(c2omegamc*(exp(c2c))/((1/(exp(c1rer))/(exp(c1rpcpd)))^((1+c2rhoc)/c2rhoc))); 
c1y     = log((exp(c1cd)+1/c1relc2*(exp(c2mc)))/(1-c1sharegy*exp(c1shockgy))); 
c1g     = log(c1sharegy*exp(c1y)*exp(c1shockgy)); 
c1hp    = log(((1+(c1nup))/(c1nup) ... 
            *c1chi0*(exp(c1y)/exp(c1z))^(c1chi)/((exp(c1c))^(-c1sigma)/exp(c1rpcpd))*exp(c1y)) ... 
            /(1- c1xip*c1beta*(((exp(c1c))^(-c1sigma)/exp(c1rpcpd)))/(((exp(c1c))^(-c1sigma)/exp(c1rpcpd)))*((c1pistar)/(c1pid))^(-(1 ... 
            +(c1nup))/(c1nup)*(1+c1chi)))); 
c1gp    = log(((1+c1taup*exp(c1taupvar))/(c1nup)*(exp(c1y)))/(1-c1xip*c1beta*(((exp(c1c))^(-c1sigma)/exp(c1rpcpd)))/(((exp(c1c))^(-c1sigma)/exp(c1rpcpd)))*((c1pistar)/(c1pid))^(1 ... 
            -(1+(c1nup))/(c1nup)))); 
c1pstar = log(((exp(c1hp))/(exp(c1gp)))^(1/((1+(1+c1nup)/(c1nup)*c1chi)))); 
c2y     = log(((exp(c2cd))+(c1relc2*(exp(c1mc))))/(1-c2sharegy*(exp(c2shockgy)))); 
c2g     = log(c2sharegy*(exp(c2y))*(exp(c2shockgy))); 
c2hp    = log(( (1+(c2nup))/(c2nup) ... 
            *c2chi0*(exp(c2y)/exp(c2z))^c2chi/((exp(c2c))^(-c2sigma)/exp(c2rpcpd))*exp(c2y))/ ... 
            (1-c2xip*c2beta*(((exp(c2c))^(-c2sigma)/exp(c2rpcpd)))/(((exp(c2c))^(-c2sigma)/exp(c2rpcpd)))*((c2pistar)/(c2pid))^(-(1 ... 
            +(c2nup))/(c2nup)*(1+c2chi))));  
c2gp    = log(((1+c2taup*exp(c2taupvar))/(c2nup)*(exp(c2y)))/(1-c2xip*c2beta*(((exp(c2c))^(-c2sigma)/exp(c2rpcpd)))/(((exp(c2c))^(-c2sigma)/exp(c2rpcpd)))*((c2pistar)/(c2pid))^(1 ... 
            -(1+(c2nup))/(c2nup))));    
c2pstar  = log(( (exp(c2hp))/(exp(c2gp)))^(1/((1+(1+c2nup)/(c2nup)*c2chi)))); 

%not needed to eval fval               
c1deltad    = log(((1-c1xip)*((exp(c1pstar)))^(-(1+(c1nup))/(c1nup)*(1+c1chi)))/ ... 
                (1-c1xip*((c1pistar)/(c1pid))^(-(1+(c1nup))/(c1nup)*(1+c1chi)))); 
c2deltad    = log(((1-c2xip)*((exp(c2pstar)))^(-(1+(c2nup))/(c2nup)*(1+c2chi)))/ ... 
                   (1-c2xip*((c2pistar)/(c2pid))^(-(1+(c2nup))/(c2nup)*(1+c2chi)))); 
c1rn        = 1/( c1beta*(((exp(c1c))^(-c1sigma)/exp(c1rpcpd)))/(((exp(c1c))^(-c1sigma)/exp(c1rpcpd)))/c1pid)-1;  
c2rn        = 1/( c2beta*(((exp(c2c))^(-c2sigma)/exp(c2rpcpd)))/(((exp(c2c))^(-c2sigma)/exp(c2rpcpd)))/c2pid)-1; 
Util1       = ( (exp(c1c))^(1-c1sigma)/(1-c1sigma)-c1chi0*(exp(c1z))^(-c1chi)*(exp(c1y))^(1+c1chi)/(1+c1chi)*exp(c1deltad)); 
Util2       = ( (exp(c2c))^(1-c2sigma)/(1-c2sigma)-c2chi0*(exp(c2z))^(-c2chi)*(exp(c2y))^(1+c2chi)/(1+c2chi)*exp(c2deltad)); 
Util        = omega_welf1*Util1+omega_welf2*Util2; 
Welfcon1    = ( (exp(c1c))^(1-c1sigma)/(1-c1sigma))/(1-c1beta); 
Welf        = Util/(1-nbeta); 
Welf1       = Util1/(1-nbeta); 
Welf2       = Util2/(1-nbeta); 
c1deltadlag = c1deltad;  
c1rnlag     = c1rn;  
c1rpcpdlag  = c1rpcpd;  
c2deltadlag = c2deltad;  
c2rnlag     = c2rn;  
c2rpcpdlag  = c2rpcpd;    
c1dcore     = ( exp(c1rpcpd)/exp(c1rpcpdlag)*c1pid); 
c2dcore     = ( exp(c2rpcpd)/exp(c2rpcpdlag)*c2pid); 


%% conditions to be fullfilled in equilibrium (check consistency of current variablesguess)
fval(1) = exp(c1c) - ( ((c1omegac)^(c1rhoc/(1+c1rhoc))*exp(c1cd)^(1/(1+c1rhoc))+(c1omegamc)^(c1rhoc/(1+c1rhoc))*(exp(c1mc))^(1/(1+c1rhoc)))^(1+c1rhoc)); 
fval(2) = (exp(c2c)) - ( ((c2omegac)^(c2rhoc/(1+c2rhoc))*(exp(c2cd))^(1/(1+c2rhoc))+(c2omegamc)^(c2rhoc/(1+c2rhoc))*((exp(c2mc)))^(1/(1+c2rhoc)))^(1+c2rhoc)); 
fval(3) = (1-c1xip)*((exp(c1pstar)))^(-1/(c1nup))+c1xip*((c1pistar)/(c1pid))^(-1/(c1nup)) - ( 1); 
fval(4) = (1-c2xip)*((exp(c2pstar)))^(-1/(c2nup))+c2xip*((c2pistar)/(c2pid))^(-1/(c2nup)) - ( 1); 


%% store current values in ys_conditional

for indxxx = 1:1:size(M_.endo_names_conditional,1)
    ys_conditional(indxxx,1) = eval(M_.endo_names_conditional(indxxx,:));
end


end



