% All-bin binary beliefs (enables bin x totbin noise averaging over full dataset)
%
% Same as binary_probs.m but uses all bins instead of 3 bins around
% the resolving state
%
% Assumes in workspace: rnProb (sparse), indicators_t_T, xCdf, xProb,
%   priceDateMat, exDateMat, priceDates, exDates, priceVals,
%   set_date, set_price, zc_dates, zc_rate, zc_dtm, last_ex

%% Define variables

initialProbCutoff=0;
dontUseIfSumLessThan=.05;
dropWhenThisFarFromResolution=300;

%% Aggregate bins

xCdfPctAgg = [-10000 -.2:.05:.2 10000];
rn1=size(indicators_t_T,1);
rn3 = min(last_ex,size(indicators_t_T,2));
xCdfPctIndex = NaN(rn3,length(xCdfPctAgg));

% Collect dates with AM settlements
unique_amexdates = unique(set_date);

% Empty matrix
rnProbAgg=NaN(rn1,length(xCdfPctAgg)-1,rn3);

% For each contract, classify raw bins into percent return bins
for exDateIndex=1:rn3

    exDate = exDates(exDateIndex);

    % Date indices on which this contract traded
    priceDateIndexes=find(isfinite(priceDateMat(:,exDateIndex)));

    % Cut off if starts trading far before expiration
    if length(priceDateIndexes)>dropWhenThisFarFromResolution
        initPriceDateIndex = max(priceDateIndexes(end-dropWhenThisFarFromResolution), priceDateIndexes(1));
    else
        initPriceDateIndex = priceDateIndexes(1);
    end
    startDate = priceDates(initPriceDateIndex);

    % Find current price of contract
    startPrice = priceVals(initPriceDateIndex);

    % AM expiration?
    am_exp = length(find(unique_amexdates==exDate,1));

    % Calculate time to maturity in years
    if am_exp==1
        opt_dtm = datenum(exDate)-datenum(startDate)-1;
    else
        opt_dtm = datenum(exDate)-datenum(startDate);
    end
    opt_ytm = opt_dtm/365;

    % Calculate risk-free rate on start date
    temp_zc_date = startDate;
    while length(find(zc_dates==temp_zc_date,2,'last'))<2
        temp_zc_date = temp_zc_date-1;
    end
    select_zc_indices = zc_dates==temp_zc_date;
    zc_curve = zc_rate(select_zc_indices)/100;
    zc_curve_dtm = zc_dtm(select_zc_indices);
    zc_interp = max(interp1(zc_curve_dtm,zc_curve,opt_dtm,'linear','extrap'),0);

    % Annualized log excess return of each bin in raw data
    pctReturnsxCdf = NaN(1,length(xCdf));
    for i=1:length(xCdf)
        pctReturnsxCdf(i) = log(xCdf(i)/startPrice)-opt_ytm*zc_interp;
    end

    % Translate to indices of xCdfPctAgg
    for i=1:length(xCdfPctAgg)
        returnDiff = abs(pctReturnsxCdf-xCdfPctAgg(i));
        [tmp,pctIndex] = min(returnDiff);
        xCdfPctIndex(exDateIndex,i) = pctIndex;
    end

    % Sum up probabilities in rnProb using linear indexing for sparse matrix
    pdi2=(initPriceDateIndex:priceDateIndexes(end))';
    inds_extract = sub2ind([rn1 rn3],pdi2,exDateIndex*ones(size(pdi2)));
    loadprob = NaN(length(inds_extract),length(xProb));
    for i=1:length(inds_extract)
        induse = inds_extract(i);
        tmpload = full(rnProb(:,induse)');
        if sum(tmpload,"omitnan")>0
            tmpload(isnan(tmpload))=0;
            loadprob(i,:)=tmpload;
        end
    end

    for i=1:(length(xCdfPctAgg)-1)
        ind2 = xCdfPctIndex(exDateIndex,i):(xCdfPctIndex(exDateIndex,i+1)-1);
        rnProbAgg(pdi2,i,exDateIndex)=sum(loadprob(:,ind2,:),2);
    end

end


%% Get binary contract details (ALL bins)

binaryProbs=NaN(rn3,size(rnProbAgg,2)-1,dropWhenThisFarFromResolution+1);
twoBinTotalProbs=NaN(rn3,size(rnProbAgg,2)-1,dropWhenThisFarFromResolution+1);
binaryProbStartOfContract=NaN(rn3,size(rnProbAgg,2)-1);
binaryProbEndOfContract=NaN(rn3,size(rnProbAgg,2)-1);
binaryDates=NaT(rn3,dropWhenThisFarFromResolution+1);

    % Loop over contracts
    for exDateIndex=1:rn3

        exDateIndex;
        exDates(exDateIndex);

        priceDateIndexes=find(isfinite(priceDateMat(:,exDateIndex)));

        % Use ALL bins (not confined to resolving states)
        initPriceDateIndex = priceDateIndexes(1);
        initPriceDateIndexIndex = 1;
        binIndexes=1:size(rnProbAgg,2);

        % Extract probability data for this contract
        probData=rnProbAgg(priceDateIndexes(initPriceDateIndexIndex:end),binIndexes,exDateIndex);
        probData(probData<0)=0;
        datesToUse = priceDates(priceDateIndexes(initPriceDateIndexIndex:end));

        % Compute conditional probabilities for adjacent bin pairs
        for i=1:length(binIndexes)-1
            binaryProb=probData(:,i)./sum(probData(:,i:i+1),2);
            twoBinTotalProb = sum(probData(:,i:i+1),2);

            bothLessThan=find(probData(:,i)+probData(:,i+1)<dontUseIfSumLessThan);
            if ~isempty(bothLessThan)
                binaryProb(bothLessThan)=NaN;
            end

            if isempty(binaryProb)
                continue
            end

            % Store in reverse time order (t=1 is resolution)
            diffFinalEnd = length(priceDateIndexes(initPriceDateIndexIndex:end))-length(binaryProb);
            firstObservationIndex = 1+diffFinalEnd;
            finalObservationIndex = firstObservationIndex+length(binaryProb)-1;
            binaryProbs(exDateIndex,binIndexes(i),firstObservationIndex:finalObservationIndex)=flipud(binaryProb);
            twoBinTotalProbs(exDateIndex,binIndexes(i),firstObservationIndex:finalObservationIndex)=flipud(twoBinTotalProb);

            try
                binaryProbStartOfContract(exDateIndex,binIndexes(i))=binaryProb(find(~isnan(binaryProb),1,'first'));
            catch
                binaryProbStartOfContract(exDateIndex,binIndexes(i))=NaN;
            end
            binaryProbEndOfContract(exDateIndex,binIndexes(i))=binaryProb(end);
            binaryDates(exDateIndex,firstObservationIndex:finalObservationIndex)=flipud(datesToUse);

        end

    end
