#!/usr/bin/python3

# Stock Options Analyze - Stock market options analysis.
# Copyright (C) 2017 Gordon Irlam
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

'''
Description:

    Stock Options Analyze analyzes stock market options from an
    Effective Altruism investing perspective.

    It requires several data files, which can be downloaded from Yahoo
    Finance and downloaded or purchased from CBOE:

        ^GSPC.csv - daily S&P 500 index, Yahoo Finance.
        ^VIX.csv - daily VIX implied volatility index, Yahoo Finance.
        Optsum_2016-01-04.csv - dated ^SPX options chains, CBOE
            Datashop Optsum sample file.
        purchased/UnderlyingOptionsEODCalcs_2016-01-04.csv - dated
            ^SPX options chains and Greeks, CBOE.
        purchased/UnderlyingOptionsEODQuotes_<date>.csv - dated ^SPX
            options chains, CBOE.

    It also requires the module:

        spia.py - https://github.com/gordoni/aacalc/blob/master/web/aacalc/spia.py

    for obtaining the risk free rate, which in turn requires yield
    curve data files downloaded from the Treasury using:

        fetch_yield_curve - https://github.com/gordoni/aacalc/blob/master/web/fetch_yield_curve

    The program prints some results and outputs data files that can
    then be plotted using GnuPlot. See the code for the definition of
    each data file.

Version:

    0.1.
'''

# Incase want to use pypy.
from __future__ import print_function
from __future__ import division

from csv import reader, writer
from datetime import datetime, timedelta
from functools import reduce
from math import ceil, exp, floor, isnan, log, pi, sqrt
from operator import mul
from os.path import expanduser, normpath
from sys import path

from scipy.interpolate import interp1d
from scipy.stats import norm

p = normpath(expanduser('~/aacalc/web/aacalc'))
if p not in path:
    path.append(p)

from spia import YieldCurve

calendar_days_per_year = 365
trading_days_per_year = 252
trading_days_per_month = 21
months_per_year = 12
assert(trading_days_per_month * months_per_year == trading_days_per_year)

def mean(l):
    return sum(l) / len(l)

def var(l):
    m = mean(l)
    return sum((x - m) ** 2 for x in l) / (len(l) - 1)

def stdev(l):
    return sqrt(var(l))

def median(l):
    return sorted(l)[int((len(l) - 1) / 2)]

def geomean(l):
    return reduce(mul, l, 1) ** (1 / len(l))

def get_csv(s):

    with open(s + '.csv') as f:
        return list(reader(f))

def get_index(sym):

    index = get_csv(sym)
    head = index.pop(0)
    assert(head[0] == 'Date' and head[5] == 'Adj Close')
    index = [(date, float(s)) for (date, _, _, _, _, s, _) in index if s != 'null']

    return dict(index)

def change(index, interval):

    changes = {}
    dates = sorted(index)
    for i in range(0, len(dates) - interval):
        changes[dates[i]] = index[dates[i + interval]] / index[dates[i]]

    return changes

index_gspc = get_index('^GSPC')
daily_log_rets = {date: log(ret) for date, ret in change(index_gspc, 1).items()}

def index_volatility(date_low, date_high):

    return stdev(tuple(log_ret for (date, log_ret) in daily_log_rets.items() if date_low <= date < date_high)) * sqrt(trading_days_per_year)

def analyze_volatility():

    global daily_log_rets

    def log_vol(log_rets, period, stride):

        lv = {}
        dates = sorted(log_rets)
        for i in range(0, len(dates) - period):
            lv[dates[i]] = stdev(tuple(log_rets[dates[i + j * stride]] for j in range(int(period / stride))))

        return lv

    daily_log_vols_annual = log_vol(daily_log_rets, trading_days_per_year, 1)
    with open('stock_options-annual.csv', 'w') as f:
        w = writer(f)
        for date in sorted(daily_log_vols_annual):
            w.writerow([date, daily_log_vols_annual[date] * sqrt(trading_days_per_year)])

    daily_log_vols_monthly = log_vol(daily_log_rets, trading_days_per_month, 1)
    with open('stock_options-monthly.csv', 'w') as f:
        w = writer(f)
        for date in sorted(daily_log_vols_monthly):
            w.writerow([date, daily_log_vols_monthly[date] * sqrt(trading_days_per_year)])

    index_full_years = dict(index_gspc)

    daily_log_rets = {date: log(ret) for date, ret in change(index_full_years, 1).items()}
    monthly_log_rets = {date: log(ret) for date, ret in change(index_full_years, trading_days_per_month).items()}
    annual_log_rets = {date: log(ret) for date, ret in change(index_full_years, trading_days_per_year).items()}

    daily_log_vols = log_vol(daily_log_rets, len(daily_log_rets) - 1, 1)
    monthly_log_vols = log_vol(monthly_log_rets, len(monthly_log_rets) - trading_days_per_month, trading_days_per_month)
    annual_log_vols = log_vol(annual_log_rets, len(annual_log_rets) - trading_days_per_year, trading_days_per_year)

    # Supposedly unlike returns, volatility show no difference when measured over daily, monthly, or annual returns.
    print('daily log vol:', mean(daily_log_vols.values()) * sqrt(trading_days_per_year))
    print('monthly log vol:', mean(monthly_log_vols.values()) * sqrt(months_per_year))
    print('annual log vol:', mean(annual_log_vols.values()))

    mu = mean(annual_log_rets.values())
    sigma = mean(annual_log_vols.values()) # Used by Black-Scholes.

    mn = exp(mu + sigma ** 2 / 2)
    vl = sqrt((exp(sigma ** 2) - 1) * exp(2 * mu + sigma ** 2)) # Not used by Black-Scholes.

    print('mu, sigma:', mu, sigma)
    print('mean, vol:', mn, vl)

    index_vix = get_index('^VIX')
    with open('stock_options-vix.csv', 'w') as f:
        w = writer(f)
        for date in sorted(index_vix):
            w.writerow([date, index_vix[date] / 100])

    print('daily log vol ^SPX 2012-2016:', index_volatility('2012', '2017'))
    print('mean VIX 2012-2016:', mean(tuple(vix / 100 for (date, vix) in index_vix.items() if '2012' < date < '2017'))) # VIX overesitimates volatility by around 25%.

def analyze_call_options():

    def norm_pdf(x):

        return float(norm.pdf(x)) # De-numpify.

    def norm_cdf(x):

        return float(norm.cdf(x)) # De-numpify.

    def black_scholes(years, type, strike, rf, div, sigma, underlying_price):

        forward = underlying_price * exp((rf - div) * years)

        try:

            #d1 = (log(underlying_price / strike) + (rf - div + sigma ** 2 / 2) * years) / (sigma * sqrt(years))
            d1 = (log(forward / strike) + (sigma ** 2 / 2) * years) / (sigma * sqrt(years))

        except (OverflowError, ZeroDivisionError):

            if type == 'call':
                return max(0, (forward - strike) * exp(- rf * years))
                #return max(0, underlying_price * exp(- div * years) - strike * exp(- rf * years))
            elif type == 'put':
                return max(0, (strike - forward) * exp(- rf * years))
                #return max(0, strike * exp(- rf * years) - underlying_price * exp(- div * years))

        else:

            d2 = d1 - sigma * sqrt(years)

            if type == 'call':
                return exp(- rf * years) * (forward * norm_cdf(d1) - strike * norm_cdf(d2))
                #return underlying_price * exp(- div * years) * norm_cdf(d1) - strike * exp(- rf * years) * norm_cdf(d2)
            elif type == 'put':
                return exp(- rf * years) * (strike * norm_cdf(- d2) - forward * norm_cdf(- d1))
                #return strike * exp(- rf * years) * norm_cdf(- d2) - underlying_price * exp(- div * years) * norm_cdf(- d1)

        assert(False)

    def delta(years, type, strike, rf, div, sigma, underlying_price):

        forward = underlying_price * exp((rf - div) * years)

        try:
            #d1 = (log(underlying_price / strike) + (rf - div + sigma ** 2 / 2) * years) / (sigma * sqrt(years))
            d1 = (log(forward / strike) + (sigma ** 2 / 2) * years) / (sigma * sqrt(years))
        except (OverflowError, ZeroDivisionError):
            c_delta = 1
        else:
            c_delta = norm_cdf(d1)

        if type == 'call':
            return exp(- div * years) * c_delta
        elif type == 'put':
            return exp(- div * years) * (c_delta - 1)

        assert(False)

    def vega(years, strike, rf, div, sigma, underlying_price):

        forward = underlying_price * exp((rf - div) * years)

        try:
            #d1 = (log(underlying_price / strike) + (rf - div + sigma ** 2 / 2) * years) / (sigma * sqrt(years))
            d1 = (log(forward / strike) + sigma ** 2 / 2 * years) / (sigma * sqrt(years))
        except (OverflowError, ZeroDivisionError):
            return 0
        else:
            return underlying_price * exp(- rf * years) * norm_pdf(d1) * sqrt(years)
            #return underlying_price * exp(- div * years) * norm_pdf(d1) * sqrt(years)

        assert(False)

    class Option(object):

        def __init__(self, root, quote_date, expiration_date, type, strike, open_interest, bid, ask, underlying_bid, underlying_ask):

            self.root = root
            self.quote_date = quote_date
            self.expiration_date = expiration_date
            self.type = type
            self.strike = strike
            self.open_interest = open_interest
            self.bid = bid
            self.ask = ask
            self.underlying_bid = underlying_bid
            self.underlying_ask = underlying_ask

            if root in ('SPX', 'SZU', 'SPL', 'SZJ', 'SZT', 'SZV', 'SXJ', 'SXG', 'SYZ', 'SYF', 'SYG', 'SPZ', 'SXB', 'SPQ', 'SPT', 'SZP', 'SXY', 'SXZ', 'SXM', 'SPB', 'SPV', 'SVP', 'SYV'):
                adjust = -0.5 # AM settled. Ceratin for SPX, not sure about the rest.
            elif root in ('SPXW', 'SPXPM', 'SPXQ', 'JXD', 'SZQ', 'SAQ', 'SLQ', 'SQP', 'SKQ', 'QSE', 'QSZ', 'QZQ', 'SQG', 'SYU', 'SZD'):
                adjust = 0 # PM settled. Certain for SPXW and SPXPM not sure about the rest.
            else:
                print('Bad root:', quote_date, root)
                assert(False)
            self.years = ((datetime.strptime(expiration_date, '%Y-%m-%d') - datetime.strptime(quote_date, '%Y-%m-%d')).days + adjust) / calendar_days_per_year

        def price(self, which = 'mid'):

            if which == 'bid':
                return self.bid
            elif which == 'ask':
                return self.ask
            elif which == 'mid':
                return (self.bid + self.ask) / 2

        def underlying_price(self, which = 'mid'):

            if which == 'bid':
                return self.underlying_bid
            elif which == 'ask':
                return self.underlying_ask
            elif which == 'mid':
                return (self.underlying_bid + self.underlying_ask) / 2

        def compute_iv(self, rf, div, which = 'mid'):

            price = self.price(which)

            if price == 0:
                return 0

            sigma = sigma_init = 0.15 # Initial guess,

            sigma_lo = 0
            sigma_hi = sigma_hi_init = 100

            max_count = 50
            for count in range(max_count):

                bs = black_scholes(self.years, self.type, self.strike, rf, div, sigma, self.underlying_price())
                v = vega(self.years, self.strike, rf, div, sigma, self.underlying_price())

                if abs(bs - price) < 1e-12 * price:
                    break

                if bs > price:
                    sigma_hi = min(sigma_hi, sigma)
                else:
                    sigma_lo = max(sigma_lo, sigma)

                try:
                    sigma -= (bs - price) / v
                except ZeroDivisionError:
                    if bs > price:
                        if sigma == 0:
                            sigma = float('-inf')
                            break
                        else:
                            sigma = 0
                    else:
                        sigma = float('inf')
                        break
                if sigma == float('inf'):
                    break

                sigma = max(0, sigma)

                if not sigma_lo <= sigma <= sigma_hi:
                    sigma = (sigma_lo + sigma_hi) / 2

            if sigma != float('inf') and sigma_hi_init - sigma < 1e-12 * sigma_hi_init:
                sigma = float('inf')

            assert(count < max_count - 1)

            return sigma

    class Strike(object):

        def __init__(self, strike):

            self.strike = strike

        def compute_idiv(self, rf, which = 'mid'):

            assert(self.call.years == self.put.years)
            assert(self.call.underlying_price() == self.put.underlying_price())
            try:
                idiv = - log((self.call.price(which) - self.put.price(which) + self.strike * exp(- rf * self.call.years)) / self.call.underlying_price()) / self.call.years
            except ZeroDivisionError:
                idiv = 0

                # Implied dividend by re-arranging put-call parity.
            #assert(idiv >= 0)
                # Fails for SPX 2016-01-15 expiration date options, expire very soon after 2016-01-04 quote date.
                #
                # The call at 500 strike is 1515.00 bid and 1517.30 ask when underlying in 2012.66 at close. Bid is too high. Don't understand why.
                #     Closing calls show increasing implied dividends as a function of contract duration.
                #     Option price is always postive.
                #
                # The call at 500 strike is 1499.50 bid and 1503.70 ask when underlying in 2003.32 at 15:45. This is OK.
                #     Closing calls show decreasing implied dividends as a function of contract duration.
                #     Option price premium becomes negative for deep in-the-money calls.
                #
                # If calls were overpriced you would want to short a call, buy a put, buy the stock, receive the discounted strike, and receive stock loaning fees.
                # If calls were underpriced you would want to buying a call, shorting a put, short the stock, pay the discounted strike, and pay stock borrowing fees.
                #
                # Implied dividends are a fiction that ensure put-call parity, and thus the volatility surfaces for puts and calls are equal.
                # Implied dividends incorporate information about the mispricing of options from the underlying price.
                #
                # Put-call parity would be violated by hard to borrow stocks.
                # The price paid to borrow hard to borrow stocks in order to short them and maintain put-call parity is effectively a negative dividend.
                # In the case of close to the expiration options there might not be time to receive any stock loaning fees.
                # But the assumption that you will receive stock loan fees is baked into the price of the stock,
                # meaning that calls would remain overpriced and arbitraige would fail.
                # Vanguard's S&P 500 fund only recevied $10m from security lending on $280b in assets, so barking up the wrong tree.
                #
                # Perhaps brokerage commissions need to be taken into account.
                # CBOE charges around $0.90 per contract for the market taker and refunds around $0.80 per contract for the market maker for trades on C2.
                # Consider a fee of $1 per contract.
                # A contract is 100 times the index, so this is 5e-6 of the underlying price,
                # or 1e-5 of the underlying price for a short term, deep in the money call at half the underlying price.
                # If this cost has to be recovered in 3 months it would be a 4e-5 drain, or 0.004%.
                # So this doesn't explain it either.

            return idiv

    class Chain(object):

        def __init__(self, quote_date):

            self.quote_date = quote_date
            self.strikes = {}

    def analyze_date(date, div_expected, src_prefix = 'purchased/UnderlyingOptionsEODQuotes_'):

        durations = {'3 month': 0.25, '6 month': 0.5, '1 year': 1.0}

        type_map = {'C': 'call', 'c': 'call', 'P': 'put', 'p': 'put'}

        options = {}
        nearest = {}
        distance = dict((desc, float('inf')) for desc in durations.keys())
        fname = src_prefix + date + '.csv'
        try:
            fh = open(fname)
        except IOError:
            print('Skipping non-existent file:', fname)
            return {}
        with fh as f:
            csv = reader(f)
            head = next(csv)
            fields = dict((field, i) for i, field in enumerate(head))
            for opt in csv:
                if opt[fields['underlying_symbol']] == '^SPX' and not opt[fields['root']] in ('BSF', 'BSK', 'BSZ', 'VSTRP', 'SRO', 'SPXW', 'SPXPM'):
                    # BSF, BSK, BSZ, VSTRP, and SRO are wierd option types.
                    # SPXW and probably the related historical SPXPM have wide and volatile spreads for 3 month and longer dates.
                    if src_prefix == 'Optsum_':
                        option = Option(opt[fields['root']], opt[fields['quote_date']], opt[fields['expiry']], type_map[opt[fields['type']]], float(opt[fields['strike']]), int(opt[fields['open_interest']]), float(opt[fields['last_bid_price']]), float(opt[fields['last_ask_price']]), float(opt[fields['underlying_close']]), float(opt[fields['underlying_close']]))
                    else:
                        #when = 'eod'
                        when = '1545'
                        option = Option(opt[fields['root']], opt[fields['quote_date']], opt[fields['expiration']], type_map[opt[fields['option_type']]], float(opt[fields['strike']]), int(opt[fields['open_interest']]), float(opt[fields['bid_' + when]]), float(opt[fields['ask_' + when]]), float(opt[fields['underlying_bid_' + when]]), float(opt[fields['underlying_ask_' + when]]))
                            # Don't understand how CBOE computes implied_underlying_price_1545, a forward price that is constant for a given chain.
                            # implied_underlying_price_1545 is not set to maintain put-call parity, so CBOE has different call and put implied volatilities.
                            # implied_underlying_price_1545 / underlying_bid_1545 is not an exponentially decreasing function of duration (i.e. they have variable dividends).
                            # Even given implied_underlying_price_1545, don't understand how CBOE computes implied_volatility_1545.
                    if option.ask == 0:
                        continue
                    expiration_date = option.expiration_date
                    for desc, duration in durations.items():
                        dist = abs(option.years - duration)
                        if (dist, expiration_date) < (distance[desc], nearest.get(desc)):
                            nearest[desc] = expiration_date
                            distance[desc] = dist
                    try:
                        chain = options[expiration_date]
                    except KeyError:
                        chain = options[expiration_date] = Chain(option.quote_date)
                    try:
                        strike = chain.strikes[option.strike]
                    except KeyError:
                        strike = chain.strikes[option.strike] = Strike(option.strike)
                    try:
                        existing_option = getattr(strike, option.type)
                        if (option.quote_date, option.expiration_date, option.strike) == ('2008-12-22', '2009-12-19', 725.0):
                            # Known data error. Second data is invalid.
                            continue
                        raise Exception('Multiple options with same quote date, expiration date, and strike: ' + option.quote_date + ' ' + option.expiration_date + ' ' + str(option.strike) + ' ' + existing_option.root + ' ' + option.root)
                    except AttributeError:
                        pass
                    setattr(strike, option.type, option)

        for label, dist in distance.items():
            if dist > 1 / months_per_year / 2:
                print('Distance of chain from target duration is large:', date, label, dist * calendar_days_per_year)

        labels = dict((expiration, label) for label, expiration in nearest.items())

        for expiration in tuple(options.keys()):
            if expiration not in labels:
                del options[expiration]

        atm_iv = {}
        with open('stock_options-implied_volatility-' + date + '.csv', 'w') as f:
            csv = writer(f)
            quote_date = tuple(options.values())[0].quote_date
            yield_curve = YieldCurve('nominal', quote_date)
            rf = log(1 + yield_curve.risk_free_rate) # Use 1 month T-bill rate as risk free rate.
            for expiration_date in sorted(options):
                label = labels[expiration_date]

                above_iv = None
                chain = options[expiration_date]
                for strike_price in sorted(chain.strikes):
                    strike = chain.strikes[strike_price]
                    try:
                        call = strike.call
                        put = strike.put
                    except AttributeError:
                        continue # Partner ask was zero.
                    assert(call.quote_date == put.quote_date)
                    assert(call.root == put.root)
                    assert(call.years == put.years)
                    assert(call.underlying_price() == put.underlying_price())
                    root = call.root
                    underlying_price = call.underlying_price()

                    idiv = strike.compute_idiv(rf)
                    strike.idiv = idiv
                    c_iv = call.compute_iv(rf, idiv)
                    p_iv = put.compute_iv(rf, idiv)
                    avg_iv = mean((c_iv, p_iv))
                    assert((abs(c_iv - p_iv) <= 1e-6 * avg_iv) or (c_iv == p_iv == float('inf') or isnan(c_iv) or isnan(p_iv) or call.years == 0))
                        # Should be equal due to put-call parity being enforced by selection of implied dividend.
                    strike.iv = avg_iv
                    csv.writerow([expiration_date, root, label, strike_price, idiv, avg_iv])

                    if strike_price <= underlying_price:
                        below = strike_price
                        below_iv = avg_iv
                    if strike_price > underlying_price and above_iv == None:
                        above = strike_price
                        above_iv = avg_iv

                atm_iv[expiration_date] = ((above - underlying_price) * below_iv + (underlying_price - below) * above_iv) / (above - below)

        results = {}
        with open('stock_options-leverage-' + date + '.csv', 'w') as f:
            csv = writer(f)
            quote_date = tuple(options.values())[0].quote_date
            vol_range_low = (datetime.strptime(quote_date, '%Y-%m-%d') - timedelta(days = calendar_days_per_year)).strftime('%Y-%m-%d')
            vol_range_high = quote_date
            vol_expected = index_volatility(vol_range_low, vol_range_high)
            yield_curve = YieldCurve('nominal', quote_date)
            rf = log(1 + yield_curve.risk_free_rate)
            for expiration_date in sorted(options):
                label = labels[expiration_date]
                atm_iv_surface = atm_iv[expiration_date]
                chain = options[expiration_date]
                overheads = []
                for strike_price in sorted(chain.strikes):
                    strike = chain.strikes[strike_price]
                    try:
                        call = strike.call
                        put = strike.put
                    except AttributeError:
                        continue # Partner ask was zero.
                    idiv_surface = strike.idiv
                    iv_surface = strike.iv
                    root = call.root
                    years = call.years
                    underlying_price = call.underlying_price()
                    try:
                        call_spread_expected = (call.ask - call.bid) / ((call.bid + call.ask) / 2)
                    except ZeroDivisionError:
                        call_expected_spread = 0
                    try:
                        put_spread_expected = (put.ask - put.bid) / ((put.bid + put.ask) / 2)
                    except ZeroDivisionError:
                        put_expected_spread = 0
                    strike_fraction = strike_price / underlying_price

                    div = div_expected
                    sigma = vol_expected
                    fair_iv = iv_surface / atm_iv_surface * sigma
                    d = delta(years, 'call', strike_fraction, rf, idiv_surface, iv_surface, 1)
                    expected_price = black_scholes(years, 'call', strike_fraction, rf, idiv_surface, iv_surface, 1)
                    fair_price = black_scholes(years, 'call', strike_fraction, rf, div, fair_iv, 1)
                    leverage = d * 1 / expected_price
                    premium = expected_price / fair_price - 1
                    annual_overhead = (premium + call_spread_expected / 2) / years
                    csv.writerow([expiration_date, root, label, strike_fraction, iv_surface, call_spread_expected, put_spread_expected, leverage, premium, annual_overhead])
                    overheads.append((leverage, annual_overhead, strike_fraction))
                overheads.sort()
                results[(label, quote_date)] = overheads

        return results

    # Next two results just intended as a usage example. Can ignore distance of chain from target warnings.
    discard = {}
    discard.update(analyze_date('2016-01-04', 0.020, 'Optsum_')) # Uses market close quotes, which have less liquidity.
    discard.update(analyze_date('2016-01-04', 0.020, 'purchased/UnderlyingOptionsEODCalcs_')) # Overwrites "Optsum_" data with 3:45pm quotes if file exists.
        # Dividend yield is 2015 year end TTM.

    results = {}
    results.update(analyze_date('2004-12-20', 0.016)) # Expected future dividend yield are all equal to TTM value as of the indicated year end.
    results.update(analyze_date('2005-12-19', 0.017))
    results.update(analyze_date('2006-12-18', 0.017))
    discard.update(analyze_date('2007-12-24', 0.018))
    discard.update(analyze_date('2008-12-22', 0.032))
    discard.update(analyze_date('2009-12-21', 0.023))
    results.update(analyze_date('2010-12-20', 0.018))
    results.update(analyze_date('2011-12-19', 0.020))
    results.update(analyze_date('2012-12-24', 0.020))
    results.update(analyze_date('2013-12-23', 0.018))
    results.update(analyze_date('2014-12-22', 0.018))
    results.update(analyze_date('2015-12-21', 0.020))
    results.update(analyze_date('2016-12-19', 0.020))

    with open('stock_options-overhead-mean.csv', 'w') as f:
        csv = writer(f)
        label_quote_dates = iter(sorted(results))
        try:
            label_quote_date = next(label_quote_dates)
        except StopIteration:
            label_quote_date = (None, None)
        prev_label, prev_quote_date = label_quote_date
        overhead_interps = []
        strike_fraction_interps = []
        while True:
            label, quote_date = label_quote_date
            if label != prev_label:
                for i in range(1000):
                    leverage = i / 100
                    overhead = 0
                    strike_fraction = 0
                    count = 0
                    for overhead_interp, strike_fraction_interp in zip(overhead_interps, strike_fraction_interps):
                        try:
                            overhead += overhead_interp(leverage)
                        except ValueError:
                            pass
                        else:
                            strike_fraction += strike_fraction_interp(leverage)
                            count += 1
                    try:
                        overhead /= count
                        strike_fraction /= count
                    except ZeroDivisionError:
                        pass
                    else:
                        csv.writerow([prev_label, leverage, overhead, strike_fraction])
                if label == None:
                    break
                prev_label = label
                overhead_interps = []
                strike_fraction_interps = []
            leverages, overheads, strike_fractions = zip(*results[label_quote_date])
            overhead_interps.append(interp1d(leverages, overheads))
            strike_fraction_interps.append(interp1d(leverages, strike_fractions))
            try:
                label_quote_date = next(label_quote_dates)
            except StopIteration:
                label_quote_date = (None, None)

analyze_volatility()
analyze_call_options()