function [BinsX,BinsY,BinsFitY]=binscatterplotOUTLIERcorrect(YY,XX,ZZ,nbins,OrderPol,ShallIplot,NamesYX)

if nargin<6
    ShallIplot=0;
end

CUToffPerc = .99;

%% Orthogonalize ZZ
Const = ones(size(YY,1),1);
ZZ = [Const ZZ];

if size(ZZ,2)>1
    [betasY]=OutlierAdjustedOLS(YY,ZZ,CUToffPerc);
    OrthY = YY - ZZ(:,2:end)*betasY(2:end);
    
    
    [betasX]=OutlierAdjustedOLS(XX,ZZ,CUToffPerc);
    OrthX = XX - ZZ(:,2:end)*betasX(2:end);
else
    OrthY = YY;
    OrthX = XX;
end

%% 
% OrthY = OrthY(:);
% OrthX = OrthX(:);

%% Take linear regression YY on XX
OrthX_all = [ones(size(OrthX))];
for pol = 1:OrderPol; 
    OrthX_all = [OrthX_all OrthX.^pol];
end

[beta,Youtadj,Xoutadj]=OutlierAdjustedOLS(OrthY,OrthX_all,CUToffPerc);
FitOrthY = Xoutadj*beta;
Xoutadj = Xoutadj(:,2);

%% Calculate Bins
[~,Sorter] = sort(Xoutadj);
OrthY = Youtadj(Sorter,1);
OrthX = Xoutadj(Sorter,:);
FitOrthY = FitOrthY(Sorter,1);

ObsPerBin = ceil(size(OrthY,1)/nbins);

for jj = 1:nbins
    if ObsPerBin*jj<=size(OrthX,1)
    BinsX(jj,1) = mean(OrthX(ObsPerBin*(jj-1)+1:ObsPerBin*jj,1));
    BinsY(jj,1) = mean(OrthY(ObsPerBin*(jj-1)+1:ObsPerBin*jj));
    BinsFitY(jj,1) = mean(FitOrthY(ObsPerBin*(jj-1)+1:ObsPerBin*jj));
    end
end


% GridX = [min(XX):(max(XX)-min(XX))/(nbins):max(XX)]';
% 
% for jj = 2:size(GridX,1)
%     Select = (XX>=GridX(jj-1) & XX<=GridX(jj));
%     BinsX(jj-1,1) = mean(OrthX(Select,2));
%     BinsY(jj-1,1) = mean(OrthY(Select));
%     BinsFitY(jj-1,1) = mean(FitOrthY(Select));
%     
% end

%% Plot?
if ShallIplot ==1
    figure
%     scatter(BinsX,BinsY,'MarkerFaceColor','b','MarkerEdgeColor','b');
    plot(BinsX,BinsY,'MarkerFaceColor','b','MarkerEdgeColor','b','MarkerSize',8,'Marker','o','LineStyle','none');
    hold on; 
    plot(BinsX,BinsFitY,'color','g','linestyle','-','linewidth',2);
    xlabel(NamesYX{1}); 
    ylabel(NamesYX{2}); 
    xlim([min(BinsX)-.01 max(BinsX)+.01]); 
    grid on
end







function [betas,Youtadj,Xoutadj]=OutlierAdjustedOLS(Y,X,CUToffPerc)


betas = X\Y; 
res = Y-X*betas;
standres = res/std(res);

CUToff = norminv(CUToffPerc);
FindOutliers=abs(standres)>CUToff;

Youtadj = Y(FindOutliers==0,:); 
Xoutadj = X(FindOutliers==0,:); 
betas = Xoutadj\Youtadj; 



