function [xCdf,xProb,rnProb,priceDateMat,exDateMat,setVals,priceVals,csCount,optCountPreClean,optCountClean,cleaned10Pct,cleaned20Pct,indicators_t_T,indicators_index] = rn_probs_altsmooth(opt_prices,best_bid,best_offer,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,impl_volatility,volume)
% Modified to use smoothing spline (csaps) with fnxtr extrapolation, no GEV tails
%
% 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% for this)
tol = 1e-6; % Search tolerance
options_fmincon = optimoptions('fmincon','Display','off');
warning('off','stats:statrobustfit:IterationLimit');
warning('off','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:illConditionedMatrix');

% 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 = sparse(length(xProb),num_dates*num_exdates); % putting xProb as row to speed up later access
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));

    % Loop through all dates for which there are prices (now parallelized)
    corresp_dates = unique(opt_dates(opt_exdates==exdate));
    start_t = find(unique_dates==min(corresp_dates));
    end_t = find(full_unique_dates<=exdate,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)
    if am_exp==1
        end_t = find(full_unique_dates<exdate,1,'last');
    end
    for t=start_t:end_t

        date = full_unique_dates(t);

        % Calculate time to maturity
        if am_exp==1
            exdate_use_vec = datevec(exdate);
            exdate_use_vec(4) = 9;
            exdate_use_vec(5) = 30;
            date_vec = datevec(date);
            date_vec(4) = 16;
            opt_dtm = datenum(datetime(exdate_use_vec))-datenum(datetime(date_vec)); % Subtract trading day's worth of hours if option expires at open of expiration date
            opt_ytm = yearfrac(datetime(date_vec),datetime(exdate_use_vec));
        else
            opt_dtm = datenum(exdate)-datenum(date);
            opt_ytm = yearfrac(date,exdate);
        end
        if opt_ytm>1
            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,'last'))/100;

        % If this is expiration date, place all probability mass in bin in which option settled, record settlement value, and continue
        if t==end_t && abs(datenum(exdate)-datenum(date))<=4
            % Check if AM-settled
            if am_exp==1
                % If so, pull settlement value from appropriate input vector
                set_price = set_prices(find(set_dates==exdate,1));
            else
                % Otherwise, pull settlement value from end-of-day price vector
                set_price = und_price;
            end

            % Place probability one in bucket corresponding to settlement price
            rnCdf_slice = xCdf>=set_price; % 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;
            setVals(t,t_exp) = set_price;
            priceVals(t) = und_price;
            csCount(t,t_exp) = 1;

            continue;
        else
            % Record today's price otherwise
            priceVals(t) = und_price;
        end

        % Option features
        select_indices = opt_exdates==exdate & opt_dates==date;
        prices = opt_prices(select_indices);
        bids = best_bid(select_indices);
        offers = best_offer(select_indices);
        strikes = opt_strikes(select_indices);
        call_flags = opt_callflag(select_indices);
        vol_use = volume(select_indices);
        impl_vol_comp = impl_volatility(select_indices);

        % Extract zero-coupon curve and calculate zc yield corresponding to option maturity
        temp_zc_date = date;
        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 = 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 = date;
        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');
        fwd_use=fwd_interp;

        % Separate calls and puts and get rid of NaN prices and arbitrage violations
        c_ind = call_flags==1 & ~isnan(prices) & prices<=und_price;
        p_ind = call_flags==0 & ~isnan(prices);
        call_strikes = strikes(c_ind);
        put_strikes = strikes(p_ind);
        call_prices = prices(c_ind);
        call_bids = bids(c_ind);
        call_offers = offers(c_ind);
        put_prices = prices(p_ind);
        put_bids = bids(p_ind);
        put_offers = offers(p_ind);
        vol_call = vol_use(c_ind);
        vol_put = vol_use(p_ind);
        impl_vol_comp_call = impl_vol_comp(c_ind);
        impl_vol_comp_put = impl_vol_comp(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);
            call_bids = call_bids(c_ind2);
            call_offers = call_offers(c_ind2);
            vol_call = vol_call(c_ind2);
            impl_vol_comp_call = impl_vol_comp_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);
            put_bids = put_bids(p_ind2);
            put_offers = put_offers(p_ind2);
            vol_put = vol_put(p_ind2);
            impl_vol_comp_put = impl_vol_comp_put(p_ind2);
        end

        % Find overlapping strikes and select OTM options (Bliss & Panigirtzoglou method, no blending)
        [both_strikes,call_ind,put_ind] = intersect(call_strikes,put_strikes);
        Xlow = fwd_use;
        Xhigh = fwd_use;
        put_indices = [find(put_strikes<Xlow);put_ind(both_strikes>=Xlow & both_strikes<=Xhigh)];
        put_strikes_use = put_strikes(put_indices);
        put_prices_use = put_prices(put_indices);
        call_indices = [call_ind(both_strikes>=Xlow & both_strikes<=Xhigh);find(call_strikes>Xhigh)];
        call_strikes_use = call_strikes(call_indices);
        call_prices_use = call_prices(call_indices);

        % 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 no usable strikes
        if sum(call_indices)+sum(put_indices)==0
            continue;
        end

        % Dividend yield input now follows Figlewski (single yield for index)
        use_dyld = und_yld;

        % Strikes and volume
        strikes_all = [put_strikes_use;call_strikes_use(call_strikes_use>Xhigh)];
        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,put_strikes_use,zc_interp,opt_ytm,put_prices_use,lim_val,use_dyld,tol,{'put'});
            nanput = isnan(imp_vols_put);
            if max(nanput)==1
                clean_put = ~nanput;
                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(call_strikes_use>Xhigh)];
            end
        else
            imp_vols_put = [];
        end
        if ~isempty(call_strikes_use)
            imp_vols_call = blsimpv(und_price,call_strikes_use,zc_interp,opt_ytm,call_prices_use,lim_val,use_dyld,tol,{'call'});
            nancall = isnan(imp_vols_call);
            if max(nancall)==1
                clean_call = ~nancall;
                call_strikes_use = call_strikes_use(clean_call);
                imp_vols_call = imp_vols_call(clean_call);
                strikes_all = [put_strikes_use;call_strikes_use(call_strikes_use>Xhigh)];
            end
        else
            imp_vols_call = [];
        end

        imp_vols_all = [imp_vols_put;imp_vols_call];
        volume_all = [vol_put;vol_call];

        % Calculate blend weights for OTM put and call IVs
        [both_strikes,call_ind,put_ind] = intersect(call_strikes_use,put_strikes_use);
        Xlow = fwd_use;
        Xhigh = fwd_use;
        put_indices = [find(put_strikes_use<Xlow);put_ind(both_strikes>=Xlow & both_strikes<=Xhigh)];
        put_strikes_fin = put_strikes_use(put_indices);
        imp_vols_put_fin = imp_vols_put(put_indices);
        call_indices = [call_ind(both_strikes>=Xlow & both_strikes<=Xhigh);find(call_strikes_use>Xhigh)];
        call_strikes_fin = call_strikes_use(call_indices);
        imp_vols_call_fin = imp_vols_call(call_indices);

        strikes_all_fin = [put_strikes_fin;call_strikes_fin(call_strikes_fin>Xhigh)];
        put_weights = max(0,min(1,(Xhigh-strikes_all_fin)/(Xhigh-Xlow)));
        puts_to_pad = length(find(call_strikes_fin>Xhigh));
        imp_vols_put_padded = [imp_vols_put_fin;zeros(puts_to_pad,1)];
        calls_to_pad = length(find(put_strikes_fin<Xlow));
        imp_vols_call_padded = [zeros(calls_to_pad,1);imp_vols_call_fin];

        ivblend = put_weights.*imp_vols_put_padded + (1-put_weights).*imp_vols_call_padded;

        % Move on if fewer observations than needed to fit spline (more than 6); otherwise record # of obs used to fit spline
        strikes_all = strikes_all_fin;
        if length(strikes_all)<5
            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 smoothing spline
        lambda = 0.99;
        smooth_param = 1/(1+lambda); % need to reparameterize for matlab's csaps
        atmiv = interp1(strikes_all,ivblend,fwd_use,'linear','extrap');
        strdist_low = strikes_all(2)-strikes_all(1);
        strike_add_low = max(strikes_all(1)-3*strdist_low,min_indval);
        strdist_high = strikes_all(end)-strikes_all(end-1);
        strike_add_high = min(strikes_all(end)+3*strdist_high,max_indval);
        strikes_all = [strike_add_low;strikes_all;strike_add_high];
        ivblend = [ivblend(1);ivblend;ivblend(end)];
        deltafit = blsdelta(und_price,strikes_all,zc_interp,opt_ytm,atmiv,use_dyld);
        vegafit = blsvega(und_price,strikes_all,zc_interp,opt_ytm,ivblend,use_dyld);
        grad_finer = 0.1;
        delta_finer = grad_finer*2;
        evalmin = min_indval - grad_finer;
        evalmax = max_indval + grad_finer;
        x_eval_int = (evalmin:grad_finer:evalmax)';
        deltafit_full = blsdelta(und_price,x_eval_int,zc_interp,opt_ytm,atmiv,use_dyld);
        fitted_ivs_pp = csaps(deltafit,ivblend,smooth_param,[],vegafit);
        fitted_ivs_pp = fnxtr(fitted_ivs_pp);
        fitted_ivs_full = fnval(fitted_ivs_pp,deltafit_full);
        if max(isnan(fitted_ivs_full))==1
            continue;
        end

        % Convert fitted IVs back to call prices
        fitted_ivs_full(fitted_ivs_full<=0)=1e-4;
        [fitted_call,~] = blsprice(und_price,x_eval_int,zc_interp,opt_ytm,fitted_ivs_full,use_dyld);
        rnCdf_slice_int = 1+exp(zc_interp*opt_ytm)*(fitted_call(3:end)-fitted_call(1:end-2))/delta_finer;
        xCdf_int = x_eval_int(2:end-1);

        % Interpolate CDF to standard grid
        rnCdf_use = interp1(xCdf_int,rnCdf_slice_int,xCdf);
        rnPdf_slice = rnCdf_use(2:end)-rnCdf_use(1:end-1);

        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;

    end

    % Intermediate save
     if mod(t_exp,20)==0||t_exp==1
         tempname = 'Output/rnProbsAltSmoothTemp.mat';
         tempstr = sprintf(tempname);
         save(tempstr,'xCdf','xProb','rnProb','priceDateMat','exDateMat','setVals','priceVals','indicators_t_T','indicators_index','-v7.3');
     end

end


end
