% Construct intraday risk-neutral probabilities from option prices
%
% Loads intraday options data, constructs risk-neutral probabilities via rn_probs_intraday.m

%% Data import from WRDS OptionMetrics download

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

% Import if new
if new==1
    
    cd Input

    % Import csv files using Matlab auto-generated functions
    addpath('Intraday');
    d = dir('Intraday/*.csv');
    n = length(d);
    if n==0 % pseudo-data switch
        addpath('Pseudo/Intraday');
        d = dir('Pseudo/Intraday/*.csv');
        n = length(d);
    end
    num_rows = 40000000;
    date1 = NaT(num_rows,1);
    exdate = NaT(num_rows,1);
    strike_price = NaN(num_rows,1);
    cp_flag = categorical(NaN(num_rows,1));
    close_price = NaN(num_rows,1);
    volume = NaN(num_rows,1);
    best_bid = NaN(num_rows,1);
    best_offer = NaN(num_rows,1);
    underlying_bid = NaN(num_rows,1);
    underlying_ask = NaN(num_rows,1);
    row_start = 1;
    for i=1:n
        d(i).name
        [~,date1t,~,exdatet,strike_pricet,cp_flagt,~,~,~,close_pricet,  ...
          volumet,~,best_bidt,~,best_offert,underlying_bidt,            ...
           underlying_askt,~,~,~,~,~,~] = importfile7(d(i).name); % Option prices
        row_end = row_start+length(date1t)-1;
        date1(row_start:row_end) = date1t;
        exdate(row_start:row_end) = exdatet;
        strike_price(row_start:row_end) = strike_pricet;
        cp_flag(row_start:row_end) = cp_flagt;
        close_price(row_start:row_end) = close_pricet;
        volume(row_start:row_end) = volumet;
        best_bid(row_start:row_end) = best_bidt;
        best_offer(row_start:row_end) = best_offert;
        underlying_bid(row_start:row_end) = underlying_bidt;
        underlying_ask(row_start:row_end) = underlying_askt;
        clear date1t exdatet strike_pricet cp_flagt close_pricet volumet best_bidt best_offert underlying_bidt underlying_askt;
        row_start = row_end+1;
    end
    row_end
    num_rows
    date1 = date1(1:row_end);
    exdate = exdate(1:row_end);
    strike_price = strike_price(1:row_end);
    cp_flag = cp_flag(1:row_end);
    close_price = close_price(1:row_end);
    volume = volume(1:row_end);
    best_bid = best_bid(1:row_end);
    best_offer = best_offer(1:row_end);
    underlying_bid = underlying_bid(1:row_end);
    underlying_ask = underlying_ask(1:row_end);
    clear i n d num_rows row_start row_end
    date2 = date1;
    ind_price = (underlying_bid+underlying_ask)/2;
    save intraday_option_prices.mat;
    clear;
    [secid3,date3,ticker3,rate] = importfile3('div_yields.csv'); % Index dividend yield
    save div_yields.mat;
    clear;
    [date4,dtm,zc_rate] = importfile4('zero_coupons.csv'); % Zero-coupon yields
    save zero_coupons.mat;
    clear;
    [secid5,date5,dtm2,forward_price,cp_flag5,index_flag5] =            ...
                                                   importfile5('atmf.csv'); % Forward prices
    save atmf.mat;
    clear;
    [set_date,set_price] = importfile6('settlement_vals.csv'); % AM settlement prices
    save set_prices.mat;
    clear;

    cd ..
    
end
    
% Use saved data
load div_yields.mat;
load zero_coupons.mat;
load atmf.mat;
load set_prices.mat;
load intraday_option_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';
hrs = hour(date1);
mins = minute(date1);
hrs_comp = hrs+mins/60;

pos_bid = best_bid>0 & best_offer>0 & hrs_comp<=16; % throw out data after 4 pm
sortmat = [exdate_vals date_vals call_flag strike_price best_bid volume]; % Sort puts before calls
[~,ind] = sortrows(sortmat,[1 2 3 4 -5 -6]); % 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);
exdate = exdate(ind);
strike_price = strike_price(ind);
volume = volume(ind);
date_vals = date_vals(ind);
exdate_vals = exdate_vals(ind);
date2 = date2(ind);
ind_price = ind_price(ind);
[date2,ind_temp,~] = unique(date2,'last');
ind_price = ind_price(ind_temp);

%% 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_intraday(opt_prices,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,volume);
exDates = min(exDateMat,[],1)';
priceDates = min(priceDateMat,[],2);