% Compute intraday binary beliefs and estimate noise (altsmooth specification)
%
% Assumes in workspace: rnProb (sparse), indicators_t_T

addpath(genpath(pwd));

% Index bins (from daily binary beliefs produced by mat_08run_all_matlab_altsmooth)
% Selective load to avoid overwriting intraday variables already in workspace
load Output/rnProbsBinaryAltSmooth.mat exDates xCdfPctIndex xCdfPctAgg;
exDatesFull = exDates;
binIndFull = xCdfPctIndex;

% Restore intraday exDates (overwritten by rnProbsBinary load above)
load Output/rnProbsIntradayAltSmooth.mat exDates priceDates priceDateMat priceVals;
last_ex_ind = last_ex;

% Get dimensions for linear indexing of sparse rnProb
num_dates = size(indicators_t_T,1);
num_exdates = size(indicators_t_T,2);
sum_cutoff = 0.05;
exDateInds = NaN(size(exDates));
exDateIndsOwn = NaN(size(exDates));
for i=1:length(exDateInds)
    if isnat(exDates(i))
        continue;
    end
    exInd = find(exDatesFull==exDates(i));
    if isempty(exInd)
        continue;
    end
    if exInd > last_ex_ind
        continue;
    end
    exDateInds(i) = exInd;
    exDateIndsOwn(i) = i;
end

priceDays = datetime(datevec(priceDates).*[1 1 1 0 0 0]);
priceDaysUse = unique(priceDays);
exDateIndsUse = exDateInds(~isnan(exDateIndsOwn));
exDateIndsOwnUse = exDateIndsOwn(~isnan(exDateIndsOwn));
exDatesOwn = exDates(exDateIndsOwnUse);
rnProbIntraday = NaN(30,length(exDateIndsOwnUse),390,length(xCdfPctAgg)-1);
sumIntraday = NaN(30,length(exDateIndsOwnUse),390,length(xCdfPctAgg)-2);
binaryIntraday = sumIntraday;

for i=1:length(exDateIndsOwnUse)

    exDateInd = exDateIndsOwnUse(i);
    exDateIndFullData = exDateIndsUse(i);
    binIndsUse = binIndFull(exDateIndFullData,:);
    priceDateIndexes=find(isfinite(priceDateMat(:,exDateInd)));
    dayVals = unique(priceDays(priceDateIndexes));
    dayValInds = NaN(size(dayVals));
    priceDateIndsByDay = reshape(priceDateIndexes,390,length(dayVals));
    for j=1:length(dayVals)
        dayValInds(j) = find(priceDaysUse==dayVals(j));
    end
    
    ind = 0;
    for j=dayValInds'
        
        ind = ind+1;
        priceDateIndsUse = priceDateIndsByDay(:,ind);

        % Convert to linear indices for sparse rnProb access
        linInds = sub2ind([num_dates num_exdates], priceDateIndsUse, exDateInd*ones(size(priceDateIndsUse)));

        % Extract all bins for these time points (rnProb is now xProb x linearIndex)
        rnProbSlice = full(rnProb(:,linInds))'; % now size 390 x length(xProb)

        for k=1:length(xCdfPctAgg)-1

            binIndLow = binIndsUse(k);
            binIndHigh = binIndsUse(k+1)-1;
            rnProbIntraday(j,i,:,k) = sum(rnProbSlice(:,binIndLow:binIndHigh),2);
            if k>1
                sumIntraday(j,i,:,k-1) = sum(rnProbIntraday(j,i,:,k-1:k),4);
                binaryIntraday(j,i,:,k-1) = rnProbIntraday(j,i,:,k-1)./sumIntraday(j,i,:,k-1);
            end
            
        end
                
    end
    
end

binaryIntraday(sumIntraday<sum_cutoff) = NaN;
binaryIntraday(binaryIntraday>1) = 1;
binaryIntraday(binaryIntraday<0) = 0;

spec.max = 50;
jmax = 60;
jvec = 0:jmax;
spec.tol = 0.05;
spec.errwin = 3;

noiseEsts = NaN(30,length(exDateIndsOwnUse),length(xCdfPctAgg)-2,length(jvec));
altNoiseEsts = NaN(30,length(exDateIndsOwnUse),length(xCdfPctAgg)-2);
noiseEstVar = altNoiseEsts;

knUse = [2 6 6 7 4 7 5 2 27]'; % output of first-stage estimation with max 200 in ChooseKn.m
day = [];
exday = [];
bin = [];

for i=1:30
    
    for j=1:length(exDateIndsOwnUse)
        
        exDateInd = exDateIndsOwnUse(j);
        
        for k=1:length(xCdfPctAgg)-2
           
            binaryUse = squeeze(binaryIntraday(i,j,:,k));
            binaryUse = binaryUse(~isnan(binaryUse));
            if length(binaryUse)<jmax
                continue
            end
            
             day=[day i];
             exday = [exday j];
             bin=[bin k];

            kn = knUse(k);
            noiseEsts(i,j,k,:) = ReMeDI2(binaryUse,jvec,kn);
            altNoiseEsts(i,j,k) = mean(diff(binaryUse).^2)/2;
            noiseEstVar(i,j,k) = ReMeDI2CLTVarR(binaryUse,0,kn,kn,1);
            
        end
        
    end
    
end

noiseEstMeans = squeeze(nanmean(noiseEsts,[1 2]));
noiseEstStd = squeeze(nanstd(noiseEsts,1,[1 2]));
altNoiseEstMeans = squeeze(nanmean(altNoiseEsts,[1 2]));
noiseEstVarMeans = squeeze(nanmean(noiseEstVar,[1 2]));

noise1 = noiseEstMeans(:,1)';
noise2 = noiseEstMeans(:,2)';
noisestd = noiseEstStd(:,1)';
noisestdclt = (noiseEstVarMeans.^(1/2))';

noisealt = altNoiseEstMeans';
noisealt2 = noisealt+noise2;

noiseEstUse = noiseEsts(:,:,2:end-1,:);
plotmean = squeeze(nanmean(noiseEstUse,[1 2 3]));
noisemeans = plotmean;

% Save is handled by the calling script (mat_09run_all_matlab_intraday_altsmooth.m)