% Construct risk-neutral probabilities from option prices (altsmooth specification)
%
% Loads options data, constructs risk-neutral probabilities via rn_probs_altsmooth.m

%% Data import from WRDS OptionMetrics download

% Flag whether to import or use previous data
new = 0;
reload_opt = 0;

% Import if new
if new==1
    
    % Import csv files using Matlab auto-generated functions
    if reload_opt==1
        [secid,date1,exdate,cp_flag,strike_price,best_bid,best_offer,volume,...
         open_interest,impl_volatility,delta,optionid,index_flag,issuer,    ...
         exercise_style] = importfile('option_prices.csv'); % Option prices
        save Input/option_prices.mat;
        clear;
    end
    [secid2,date2,ticker,ind_price] = importfile2('index_prices.csv'); % Underlying index prices
    save Input/index_prices.mat;
    clear;
    [secid3,date3,ticker3,rate] = importfile3('div_yields.csv'); % Index dividend yield
    save Input/div_yields.mat;
    clear;
    [date4,dtm,zc_rate] = importfile4('zero_coupons.csv'); % Zero-coupon yields
    save Input/zero_coupons.mat;
    clear;
    [secid5,date5,dtm2,forward_price,cp_flag5,index_flag5] =            ...
                                                   importfile5('atmf.csv'); % Forward prices
    save Input/atmf.mat;
    clear;
    [set_date,set_price] = importfile6('settlement_vals.csv'); % AM settlement prices
    save Input/set_prices.mat;
    clear;
    
end

% Load saved data from Input folder
load Input/option_prices.mat;
load Input/index_prices.mat;
load Input/div_yields.mat;
load Input/zero_coupons.mat;
load Input/atmf.mat;
load Input/set_prices.mat;

%% Sort data

% Sort options by expiration date, date, call/put, and strike price
date_vals = datenum(date1); % Convert dates to numbers for sorting purposes
exdate_vals = datenum(exdate);
call_flag = cp_flag=='C';
pos_bid = best_bid>0 & best_offer>0 & year(exdate)<=2018;
zero_bid = best_bid<=0;
zero_bidoff = best_bid<=0|best_offer<=0;
prices_millions = length(pos_bid)/1000000;
prices_zeros = sum(zero_bidoff)/1000000;

sortmat = [exdate_vals date_vals call_flag strike_price volume]; % Sort puts before calls
[~,ind] = sortrows(sortmat,[1 2 3 4 -5]); % Create sort indices
ind = ind(pos_bid(ind)); % Keep only entries with non-zero bid prices
best_bid = best_bid(ind);
best_offer = best_offer(ind);
cp_flag = cp_flag(ind);
call_flag = call_flag(ind);
date1 = date1(ind);
delta = delta(ind);
exdate = exdate(ind);
impl_volatility = impl_volatility(ind);
strike_price = strike_price(ind)/1000; % OptionMetrics stores strikes in thousandths
volume = volume(ind);
date_vals = date_vals(ind);
exdate_vals = exdate_vals(ind);

%% Run risk-neutral function

% Set up inputs
opt_prices = (best_bid+best_offer)/2; % Prices will be midpoints between bid and ask
min_indval = 200; % Compute probabilities between S&P index values of 200 and 4000
max_indval = 4000;
delta = 4; % Buckets will have index value length of delta/2 (I'm using this delta notation to follow Malz, 2014)

% Run 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,date1,exdate,call_flag,strike_price,ind_price,date2,rate,date3,zc_rate,dtm,date4,forward_price,dtm2,date5,set_price,set_date,min_indval,max_indval,delta,impl_volatility,volume);
exDates = min(exDateMat,[],1)';
priceDates = min(priceDateMat,[],2);