% Aggregate RN distributions into binary beliefs for resolving states (5pp bins, altsmooth)
%
% Same structure as binary_probs.m; called from mat_08 and mat_10 (altsmooth specs)
%
% 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;
% Skip bins where total probability is too small for reliable estimates
dontUseIfSumLessThan=.05;
% Maximum distance from expiration (trading days)
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
    
    exDateIndex;
    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)); % initial date is either first date of contract or (expiration date-dropWhenThisFarFromResolution), whichever comes later
    else
        initPriceDateIndex = priceDateIndexes(1);
    end
    startDate = priceDates(initPriceDateIndex);
    pdi2=(initPriceDateIndex:priceDateIndexes(end))';

    % 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; % Subtract one if option expires at open of expiration date
    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); % Conduct linear interpolation to find zc rate corresponding to dtm of this set of options
    
    % 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
    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
% Dimensions: contract x bin pair x days from resolution
binaryProbs=NaN(size(rnProbAgg,3),size(rnProbAgg,2)-1,dropWhenThisFarFromResolution+1);
% Starting and ending binary probabilities by contract and bin pair
binaryProbStartOfContract=NaN(size(rnProbAgg,3),size(rnProbAgg,2)-1);
binaryProbEndOfContract=NaN(size(rnProbAgg,3),size(rnProbAgg,2)-1);
binaryDates=NaT(size(rnProbAgg,3),dropWhenThisFarFromResolution+1);
    
    % Loop over contracts
    for exDateIndex=1:size(rnProbAgg,3)

        exDateIndex;

        priceDateIndexes=find(isfinite(priceDateMat(:,exDateIndex)));
        if (initialProbCutoff>0) 
            % Only look at bins above the probability cutoff
            if length(priceDateIndexes)>dropWhenThisFarFromResolution
                initPriceDateIndex = max(priceDateIndexes(end-dropWhenThisFarFromResolution), priceDateIndexes(1)); % initial date is either first date of contract or (expiration date-dropWhenThisFarFromResolution), whichever comes later
                initPriceDateIndexIndex = max(length(priceDateIndexes)-dropWhenThisFarFromResolution,1);
            else
                initPriceDateIndex = priceDateIndexes(1);
                initPriceDateIndexIndex = 1;
            end
            binIndexes=find(rnProbAgg(initPriceDateIndex,:,exDateIndex)>=initialProbCutoff);
            % Keep the boundary observation when crossing the cutoff
            binIndexes=min(binIndexes):max(binIndexes);
        else
            % Use the three bins surrounding the resolution bin
            baseIndex = find(rnProbAgg(priceDateIndexes(end),:,exDateIndex)==1);
            if length(baseIndex)~=1
                continue;
            end
            initPriceDateIndex = priceDateIndexes(1);
            initPriceDateIndexIndex = 1;
            if baseIndex==1
                binIndexes = 1:2;
            elseif baseIndex==size(rnProbAgg,2)
                binIndexes = baseIndex-1:baseIndex;
            else
                binIndexes = baseIndex-1:baseIndex+1;
            end
        end

        % 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);
            
            % NaN out observations where combined probability is below cutoff
            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;
            % Store reversed in output matrix
            binaryProbs(exDateIndex,binIndexes(i),firstObservationIndex:finalObservationIndex)=flipud(binaryProb);
            
            % Store initial and final probabilities and dates
            binaryProbStartOfContract(exDateIndex,binIndexes(i))=binaryProb(find(~isnan(binaryProb),1,'first'));
            binaryProbEndOfContract(exDateIndex,binIndexes(i))=binaryProb(end);
            binaryDates(exDateIndex,firstObservationIndex:finalObservationIndex)=flipud(datesToUse);
            
        end
        
    end