function [xCdf,xProb,rnProb,priceDateMat,exDateMat,setVals,priceVals,csCount,optCountPreClean,optCountClean,cleaned10Pct,cleaned20Pct,indicators_t_T,indicators_index] = rn_probs_intraday(opt_prices,opt_dates,opt_exdates,opt_callflag,opt_strikes,ind_prices,ind_dates,ind_ylds,ind_yld_dates,zc_rates,zc_dtm,zc_dates,fwd_prices,fwd_dtm,fwd_dates,set_prices,set_dates,min_indval,max_indval,delta,volume)
% Intraday version: main specification
%
% Function to construct risk-neutral probabilities given option-price data using the methodology of Malz (2014) and averaging output from put prices and call prices
%
% INPUTS:
% opt_prices: Vector of call option prices
% opt_dates: Vector of dates corresponding to prices, sorted within each expiration date
% opt_exdates: Vector of expiration dates corresponding to prices
% opt_callflag: Vector of logical values indicating whether option is a call (1) or put (0)
% opt_strikes: Vector of strike prices corresponding to option prices
% ind_prices: Vector of prices for underlying index
% ind_dates: Vector of dates corresponding to index prices
% ind_ylds: Vector of dividend yields for underlying index
% ind_yld_dates: Vector of dates corresponding to dividend yields
% zc_rates: Vector of annualized zero-coupon rates
% zc_dtm: Vector of days to maturity corresponding to zero-coupon rates
% zc_dates: Vector of dates corresponding to zero-coupon rates
% fwd_prices: Vector of at-the-money forward prices for underlying index
% fwd_dtm: Vector of days to maturity corresponding to forward prices
% fwd_dates: Vector of dates corresponding to forward prices
% min_indval: Minimum underlying index value for which to calculate CDF and probability values
% max_indval: Maximum underlying index value for which to calculate CDF and probability values
% delta: Bin width for probability buckets for index
%
% OUTPUTS:
% xCdf: Vector of index values at which CDF is evaluated at all dates
% xProb: Vector of index values representing midpoints of bins over which probabilities are evaluated at all dates
% rnProb: Array of probabilities of price being within given bin of index values. rnProb(i,j,k) gives implied probability on date priceDates(i), for bin with midpoint xPdf(j), and for expiration date exDates(k)
% priceDateMat: Matrix of dates corresponding to values in rnProb; priceDateMat(i,k) gives date for all non-NaN values in rnProb(i,:,k)
% exDateMat: Matrix of option expiration dates corresponding to values in rnProb; exDateMat(i,k) gives expiration date for all non-NaN values in rnProb(i,:,k)
% setVals: Vector of option settlement values on expiration dates (AM if AM-settled, PM otherwise) corresponding to exDates vector (i.e., min(exDateMat,[],1)')
% priceVals: Vector of closing index values on each date corresponding to priceDates vector (i.e., min(priceDateMat,[],2)), including option settlement values on non-trading dates

% Set parameters
lim_val = 1; % Maximum Black-Scholes implied volatility (100%)
tol = 1e-6; % Search tolerance
options_fmincon = optimoptions('fmincon','Display','off');

% Number of option expiration dates and price dates
unique_exdates = unique(opt_exdates);
num_exdates = length(unique_exdates);
unique_dates = unique(opt_dates);
full_unique_dates = [unique_dates];
num_dates = length(full_unique_dates);
unique_amexdates = unique(set_dates);

% Empty arrays
xCdf = min_indval:delta/2:max_indval;
xProb = min_indval+delta/4:delta/2:max_indval-delta/4;
rnProb = spalloc(length(xProb),num_dates*num_exdates,1e8); % putting xProb as row to speed up later access (sparse matrices are stored columnwise, so want to run sequentially by column)
priceDateMat = datetime(NaN(num_dates,num_exdates),NaN,NaN); % Can't use NaT on old version of Matlab available on the cluster
exDateMat = datetime(NaN(num_dates,num_exdates),NaN,NaN);
setVals = NaN(num_dates,num_exdates);
priceVals = NaN(num_dates,1);
csCount = zeros(num_dates,num_exdates);
optCountPreClean = zeros(num_dates,num_exdates);
optCountClean = zeros(num_dates,num_exdates);
cleaned10Pct = zeros(num_dates,num_exdates,2);
cleaned20Pct = zeros(num_dates,num_exdates,2);
indicators_t_T = zeros(num_dates,num_exdates);
indicators_index = zeros(num_dates*num_exdates,1);

% Loop through all option expiration dates
for t_exp=1:num_exdates
    
    % Display progress
    t_exp
    
    % Calculate expiration and whether this option is AM-settled
    exdate = unique_exdates(t_exp);
    exdate
    am_exp = length(find(unique_amexdates==exdate,1));
    exdate_use_vec = datevec(exdate);
    if am_exp==1
        exdate_use_vec(4) = 9;
        exdate_use_vec(5) = 30;
    else
        exdate_use_vec(4) = 16;
        exdate_use_vec(5) = 0;
    end
    exdate_act = datetime(exdate_use_vec);
    
    % Loop through all dates for which there are prices (now parallelized)
    corresp_dates = unique(opt_dates(opt_exdates==exdate));
    start_t_vec = find(unique_dates==min(corresp_dates));
    start_t = start_t_vec(1);  % Extract scalar for compatibility
    end_t_vec = find(full_unique_dates<=exdate_act,1,'last'); % Will go until last trading day before expiration (even if no trading for this option on that day, as we'll place the resolution values as if they occurred on that last trading day)
    end_t = end_t_vec(1);  % Extract scalar for compatibility
    for t=start_t:end_t
        
        date = full_unique_dates(t);
        date_vec = datevec(date);
        opt_dtm = datenum(exdate_act)-datenum(date);
        opt_ytm = yearfrac(date,exdate_act);
        if opt_ytm>1
            continue
        end
        if datenum(datetime(date_vec.*[1 1 1 0 0 0]))==datenum(datetime(exdate_use_vec.*[1 1 1 0 0 0]))
            continue
        end
        
        % Check if current date has any prices corresponding to expiration date or if this is expiration date
        if isempty(find(corresp_dates==date,1)) && t~=end_t
            if t>length(unique_dates)
                priceDateMat(t,t_exp) = date;
                exDateMat(t,t_exp) = exdate;
            end
            continue;
        end        
        priceDateMat(t,t_exp) = date;
        exDateMat(t,t_exp) = exdate;
        
        % Underlying price and index yield
        und_price = ind_prices(find(ind_dates<=date,1,'last'));
        und_yld = ind_ylds(find(ind_yld_dates<date+1,1,'last'))/100;
        priceVals(t) = und_price;
                
        % Option features
        select_indices = opt_exdates==exdate & opt_dates==date;
        prices = opt_prices(select_indices);
        strikes = opt_strikes(select_indices);
        call_flags = opt_callflag(select_indices);
        vol_use = volume(select_indices);
                        
        % Extract zero-coupon curve and calculate zc yield corresponding to option maturity
        temp_zc_date = datetime(datevec(date).*[1 1 1 0 0 0]);
        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_rates(select_zc_indices)/100;
        zc_curve_dtm = zc_dtm(select_zc_indices);
        zc_interp_comp = 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
        
        % Extract forward-price curve and calculate forward price corresponding to option maturity
        temp_fwd_date = datetime(datevec(date).*[1 1 1 0 0 0]);
        while length(find(fwd_dates==temp_fwd_date,2,'last'))<2
            temp_fwd_date = temp_fwd_date-1;
        end
        select_fwd_indices = fwd_dates==temp_fwd_date;
        fwd_curve = fwd_prices(select_fwd_indices);
        fwd_curve_dtm = fwd_dtm(select_fwd_indices);
        if length(unique(fwd_curve_dtm))<length(fwd_curve_dtm)
            [fwd_curve_dtm,fwd_curve] = consolidator(fwd_curve_dtm,fwd_curve);
        end
        fwd_interp = interp1(fwd_curve_dtm,fwd_curve,opt_dtm,'linear','extrap');
        
        % Separate calls and puts and get rid of NaN prices
        c_ind = call_flags==1 & ~isnan(prices);
        p_ind = call_flags==0 & ~isnan(prices);
        call_strikes = strikes(c_ind);
        put_strikes = strikes(p_ind);
        call_prices = prices(c_ind);
        put_prices = prices(p_ind);
        vol_call = vol_use(c_ind);
        vol_put = vol_use(p_ind);
        
        % Consolidate call + put prices (delete duplicates)
        if length(unique(call_strikes))<length(call_strikes)
            [call_strikes,c_ind2,~] = unique(call_strikes,'stable');
            call_prices = call_prices(c_ind2);
            vol_call = vol_call(c_ind2);
        end
        if length(unique(put_strikes))<length(put_strikes)
            [put_strikes,p_ind2,~] = unique(put_strikes,'stable');
            put_prices = put_prices(p_ind2);
            vol_put = vol_put(p_ind2);
        end

        % Find strike at which call and put prices are closest and use this as the forward
        [both_strikes,call_ind,put_ind] = intersect(call_strikes,put_strikes);
        call_comp = call_prices(call_ind);
        put_comp = put_prices(put_ind);
        [~,fwd_use_ind] = min(abs(call_comp-put_comp));
        fwd_use = both_strikes(fwd_use_ind);
        if isempty(fwd_use)
            continue;
        end
        cp_use = call_comp(fwd_use_ind);
        pp_use = put_comp(fwd_use_ind);
        if cp_use<pp_use
            call_indices = call_strikes>=fwd_use;
            put_indices = put_strikes<fwd_use;
        else
            call_indices = call_strikes>fwd_use;
            put_indices = put_strikes<=fwd_use;
        end
            
        % Compute pre-cleaning option count
        optCountPreClean(t,t_exp) = sum(call_indices)+sum(put_indices);
        strikesPreClean = [put_strikes(put_indices);call_strikes(call_indices)];
        strikesPreCleanPctFwd = strikesPreClean/fwd_use;
        
        % Skip if <2 strikes for calculating dividend yield and interest rate
        if length(both_strikes)<2
            continue;
        end
        
        % Find implied dividend yield and interest rate using put-call parity for intersecting strikes and Theil-Sen (median) estimator
        [b2,a2] = TheilSen([both_strikes put_comp-call_comp]);
        use_d = a2 + und_price;
        use_dyld = -log(-(use_d/und_price-1))/opt_ytm;
        zc_interp = -log(b2)/opt_ytm;
        if und_price-use_d<0 || ~isreal(zc_interp)
            continue;
        end
        use_dyld = 0; % will subtract PV of dividends from underlying, so set yield to 0 for consistency
        
        % Use OTM strikes
        call_strikes_use = call_strikes(call_indices);
        call_prices_use = call_prices(call_indices);
        put_strikes_use = put_strikes(put_indices);
        put_prices_use = put_prices(put_indices);
        strikes_all = [put_strikes_use;call_strikes_use];
        vol_call = vol_call(call_indices);
        vol_put = vol_put(put_indices);
                
        % Calculate BS implied vols
        if ~isempty(put_strikes_use)
            imp_vols_put = blsimpv(und_price-use_d,put_strikes_use,zc_interp,opt_ytm,put_prices_use,lim_val,use_dyld,tol,{'put'});
            if max(isnan(imp_vols_put))==1
                clean_put = ~isnan(imp_vols_put);
                put_strikes_use = put_strikes_use(clean_put);
                imp_vols_put = imp_vols_put(clean_put);
                vol_put = vol_put(clean_put);
                strikes_all = [put_strikes_use;call_strikes_use];
            end
            if ~isempty(put_strikes_use)
                put_prices_as_calls = blsprice(und_price-use_d,put_strikes_use,zc_interp,opt_ytm,imp_vols_put,use_dyld);
            else
                put_prices_as_calls = [];
            end
        else
            imp_vols_put = [];
            put_prices_as_calls = [];
        end
        if ~isempty(call_strikes_use)
            imp_vols_call = blsimpv(und_price-use_d,call_strikes_use,zc_interp,opt_ytm,call_prices_use,lim_val,use_dyld,tol,{'call'});
        else
            imp_vols_call = [];
        end
        imp_vols_all = [imp_vols_put;imp_vols_call];
        prices_all = [put_prices_as_calls;call_prices_use];
        volume_all = [vol_put;vol_call];
        % Filters: uncomputable IVs, IV>100, call price greater than underlying
        clean_ind = ~isnan(imp_vols_all) & imp_vols_all<=1 & prices_all<=und_price;
        strikes_all = strikes_all(clean_ind);
        imp_vols_all = imp_vols_all(clean_ind);
        prices_all = prices_all(clean_ind);
        volume_all = volume_all(clean_ind);
                
        % Move on if fewer than 3 observations; otherwise record number of points used to fit spline
        if length(strikes_all)<3
            continue
        end
        optCountClean(t,t_exp)=length(strikes_all);
        strikesCleanPctFwd = strikes_all/fwd_use;
        strikesCleaned = setdiff(strikesPreCleanPctFwd,strikesCleanPctFwd);
        cleaned10Pct(t,t_exp,:) = [sum(abs(strikesCleaned-1)>0.1) length(strikesCleaned)];
        cleaned20Pct(t,t_exp,:) = [sum(abs(strikesCleaned-1)>0.2) length(strikesCleaned)];
        
        % Fit clamped cubic spline to implied vol values so that slopes at endpoints are zero
        strikes_all_pct = strikes_all/fwd_use;
        iv_spline = spline(strikes_all_pct',[0 imp_vols_all' 0]);
        
        % Calculate values (in percent of forward price) over which to evaluate spline
        grad_length = (delta/2)/fwd_use;
        min_pct_val = min_indval/fwd_use-grad_length;
        max_pct_val = max_indval/fwd_use+grad_length;         
        x_eval = min_pct_val:grad_length:max_pct_val;
        
        % Evaluate spline and convert back to call prices
        fitted_call = ppval(iv_spline,x_eval);
        fitted_call(x_eval<strikes_all_pct(1)) = imp_vols_all(1);
        fitted_call(x_eval>strikes_all_pct(end)) = imp_vols_all(end);
        fitted_call(fitted_call<0) = tol;
        strikes_extr = transpose(x_eval*fwd_use);
        [fitted_call,~] = blsprice(und_price-use_d,strikes_extr,zc_interp,opt_ytm,fitted_call',use_dyld);
        
        % Calculate CDF and PDF values
        rnCdf_slice = 1+exp(zc_interp*opt_ytm)*(fitted_call(3:end)-fitted_call(1:end-2))/delta; % CDF
        rnPdf_slice = (rnCdf_slice(2:end)-rnCdf_slice(1:end-1)); % Probability mass
        rnProbInd = sub2ind([num_dates num_exdates],t,t_exp); % linear index to put pdf slice into
        rnProb(:,rnProbInd) = rnPdf_slice;
        indicators_t_T(t,t_exp) = 1;
        indicators_index(rnProbInd) = 1;
        
        fittedValAtStrikes = ppval(iv_spline,strikes_all/fwd_use);
        [fittedValAtStrikes,~] = blsprice(und_price-use_d,strikes_all,zc_interp,opt_ytm,fittedValAtStrikes,use_dyld);
        rnCdfDenom = strikes_all(3:end)-strikes_all(1:end-2);
        rnCdfAtStrikes = 1+exp(zc_interp*opt_ytm)*(fittedValAtStrikes(3:end)-fittedValAtStrikes(1:end-2))./rnCdfDenom;
        rnPdfAtStrikes = rnCdfAtStrikes(2:end)-rnCdfAtStrikes(1:end-1);

    end
    
end


end
