@søren-pallesen

Hi, yes, among many other items I've created a rolling VWAP for backtrader which I've found useful. And as a challenge to myself I did it in four actual lines of code. If you use this you need to pay me by shipping me a box of your local delicacies (although not dried fish or lutefisk).

Investopedia has some good entries on VWAP and MVWAP as I recall.

'''
Author: B. Bradford 

MIT License

Copyright (c) B. Bradford

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

That they contact me for shipping information for the purpose of sending a 
local delicacy of their choice native to whatever region they are domiciled.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.


class VolumeWeightedAveragePrice(bt.Indicator):
    plotinfo = dict(subplot=False)

    params = (('period', 30), )

    alias = ('VWAP', 'VolumeWeightedAveragePrice',)
    lines = ('VWAP',)
    plotlines = dict(VWAP=dict(alpha=0.50, linestyle='-.', linewidth=2.0))



    def __init__(self):
        # Before super to ensure mixins (right-hand side in subclassing)
        # can see the assignment operation and operate on the line
        cumvol = bt.ind.SumN(self.data.volume, period = self.p.period)
        typprice = ((self.data.close + self.data.high + self.data.low)/3) * self.data.volume
        cumtypprice = bt.ind.SumN(typprice, period=self.p.period)
        self.lines[0] = cumtypprice / cumvol

        super(VolumeWeightedAveragePrice, self).__init__()

Here's an indicator (Haar wavelet with lifting) which, similar to the other wavelet implementations I posted, can essentially break a time series down to constituent frequency components. However because this is implemented via lifting it reduces boundary effects To be more compliant with Daniel R.'s vision, this does not use numpy and all calculations are done internally with no libraries.

This can be used as a crossover indicator for someone who wants an off the shelf solution. For the more advanced user I've also implemented basic thresholding of the wavelet coefficients (threshlist), or you can pick out a wavelet level in isolation (coeffpick).

Further, it also has an incremental flag which although very computationally expensive (an inefficiency that I may refine later on) shows the current bar value without computing from future bars. I don't know that this really makes as much of a difference in practice as an indicator, however.

This is implemented in Python 3.

To call:

examplelwt = bbhaarlift(self.data, coefflev=9, coeffpick = 3, blnincremental = True, threshlist=[1, 1, 1, 1, 0, 0, 0, 0, 0])

Where coeffpick and threshlist are mutually exclusive use one or the other. If both are provided, coeffpick takes precedence. The above picks out the third wavelet level out of nine total.

Because this uses lifting, 2**n bars are needed for calculation. If coefflev isn't provided, the routine uses the next lowest power of 2 bars. Ie. all the bars it can given your input data length.

And here's the code:

'''
Author: B. Bradford adapted from information provided by
        bearcave.com and Ian Kaplan.

MIT License

Copyright (c) B. Bradford

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
'''


import backtrader as bt


# input lengths must be a power of two for haar with lifting:


class bbhaarlift(bt.Indicator):

    plotinfo = dict(subplot=False)
    lines = ('haarlift', )
    params = (('coefflev', None),
              ('coeffpick', None),
              ('blnincremental', False),
              ('threshlist', [1]))      #2**coefflev data points required or 2**coefflev

    @staticmethod
    def next_power_of_2(x):
        return 1 if x == 0 else 2 ** (x.bit_length() -1)



    def once(self, start, end):
        self.blnforward = 1
        self.blninverse = 2

        #create dummy line to manipulate and desired 2** data length (max or selected)
        lstinputdata = self.data.array[:]                                   #make a copy of the list/array
        if self.params.coefflev is None:                                    #max available data length
            pow2padding = self.next_power_of_2(len(lstinputdata))           #input needs to be power of two data points
        else:                                                               #selected 2** from coefflev
            pow2padding = 2 ** self.params.coefflev
            if self.next_power_of_2(len(lstinputdata)) < pow2padding: raise ValueError("COEFF LEVEL %s GREATER THAN MAX DATA LENGTH %s AVAILABLE" % (self.params.coefflev, self.next_power_of_2(len(lstinputdata))))

        if self.params.coefflev is None:
            self.params.coefflev = pow2padding.bit_length() - 1

        # If specific coeffpick selected set it to one, and others to zero:
        if self.params.coeffpick is not None:
            self.params.threshlist = [1] * self.params.coeffpick + [0] * (
            self.params.coefflev - self.params.coeffpick)

        dpcnt = len(lstinputdata) - pow2padding                             #only do loop once at end of array if blnincremental not true

        #Count backwards through line to leave most recent value in place if incremental:
        while dpcnt > 0:
            dpcnt -= 1                                                      #decrement immediately as list slicing is zero-based
            inlist = self.data.array[dpcnt:dpcnt + pow2padding]

            ftresult = bbhaarlift.forwardtrans(self, vec=inlist)            #forward transform

            threshedresult = bbhaarlift.wvltthreshold(coeffary=ftresult, threshlist=self.params.threshlist)

            itresult = bbhaarlift.inversetrans(self, vec=threshedresult)    #inverse transform

            dummyaryout = lstinputdata[:dpcnt] + itresult + lstinputdata[dpcnt + pow2padding:]
            lstinputdata = dummyaryout                                      # self.lines[0].array[pow2padding:] + itresult

            #Only run once if not doing incremental:
            if not self.params.blnincremental: break
        self.lines[0].array = lstinputdata

    def wvltthreshold(coeffary, threshlist):

        print ("bbwaveletlifting:  len(coeffary).bit_length()", len(coeffary).bit_length())
        if len(coeffary).bit_length() < len(threshlist):
            raise ValueError('WARNING: THRESHOLD LIST LENGTH IS SMALLER THAN NUMBER OF COEFF LEVELS.')

        threshedresult = coeffary[:]
        for i in range(1, len(threshlist)):
            for x in range(2**i, 2**(i+1)):
                print ("bbwaveletlifting:  i, x, len(threshlist)  ", i, x, len(threshlist))
                threshedresult[x] = coeffary[x] * threshlist[i]
        return threshedresult


    def split(vec, N):
        start = 1
        end = N - 1

        while (start < end):
            #Could do this with list comprehension, but want to stick to the original java structure:
            for i in range(start, end, 2):
                tmp = vec[i]
                vec[i] = vec [i + 1]
                vec[i+1] = tmp
            start += 1
            end -= 1

        return vec

    def merge(self, vec, N):
        half = N >> 1
        start = half-1
        end = half

        while (start > 0):
            for i in range(start, end, 2):
                tmp = vec[i]
                vec[i] = vec[i+1]
                vec[i+1] = tmp
            start -= 1
            end += 1

    def forwardtrans(self, vec):
        '''
        The result of forwardTrans is a set of wavelet coefficients
        ordered by increasing frequency and an approximate average
        of the input data set in vec[0].  The coefficient bands
        follow this element in powers of two (e.g., 1, 2, 4, 8...).
        '''

        N = len(vec)
        print ("forwardtrans N: ", N)
        #set for counter:
        n = N
        print ("forwardtrans  n:  ", n)
        while n > 1:
            print ("forwardtrans split...")
            vec = bbhaarlift.split(vec, n)
            print("forwardtrans predict...")
            vec = bbhaarlift.predict (self, vec, n, self.blnforward)
            print("update split...")
            vec = bbhaarlift.update (self, vec, n, self.blnforward)
            n = n >> 1
        return vec

    def inversetrans(self, vec):
        N = len(vec)
        n = 2
        while n <= N:
            bbhaarlift.update(self, vec, n, self.blninverse)
            bbhaarlift.predict (self, vec, n, self.blninverse)
            bbhaarlift.merge(self, vec, n)
            n = n << 1
        return vec

    def predict(self, vec, N, direction):
        half = N >> 1
        cnt = 0

        for i in range (0, half, 1):
            predictval = float(vec[i])
            j = int(i+half)

            if (direction == self.blnforward):
                vec[j] = vec[j] - predictval
            elif (direction == self.blninverse):
                vec[j] = vec[j] + predictval
            else:
                debug.print("haar: predict: bad direction value")

        return vec

    def update(self, vec, N, direction):
        half = N >> 1

        for i in range (0, half, 1):
            j = int(i + half)
            updateVal = float(vec[j] / 2.0)

            if (direction == self.blnforward):
                vec[i] = vec[i] + updateVal
            elif (direction == self.blninverse):
                vec[i] = vec[i] - updateVal
            else:
                debug.print("haar update: bad direction value")

        return vec

Although this routine is great for its simple flat-file data caching, it does have an implementation problem that can shut down the connection before the downloading of data is completed with its sleep statements. Another issue is that, even for small downloads, the sleep statements hold up processing and when one is trying to do analysis of 100-200 securities/hour these sleep statements have to be removed as below.

Ideally, the callback will just terminate the connection correctly when the "finished" response is received from the api. So I've altered the code below. It's a bit hard to follow due to the callbacks but appears to work effectively. Please note the timeoutsecs argument should be specified as it is the maximum that the routine will wait for the download.

from ib.opt import ibConnection, message
from ib.ext.Contract import Contract
from time import time, strftime
from datetime import datetime
import pandas as pd


class IBDataCache(object):

    def _reset_data(self, host='127.0.0.1', port=4001, client_id=None):
        self._df = pd.DataFrame(columns=['Date', 'Open', 'High', 'Low', 'Close', 'Volume', 'OpenInterest'])
        self._s = pd.Series()

        #Initialize connection as long as it's not already connected:
        if (not hasattr(self, '_conn')) or (not self._conn.isConnected()):
            self._conn = ibConnection(host, port, client_id)
            self._conn.enableLogging()
            # Register the response callback function and type of data to be returned
            self._conn.register(self._error_handler, message.Error)
            self._conn.register(self._historical_data_handler, message.historicalData)
            self._conn.register(self._save_order_id, 'NextValidId')
            self._conn.register(self._nextValidId_handler, message.nextValidId)
            self._conn.connect()



    def _save_order_id(self, msg):
        self._next_valid_id = msg.orderId
        print('save order id',msg.orderId)

    def _nextValidId_handler(self, msg):
        print("nextValidId_handler: ", msg)
        self.inner(sec_type=self.req.m_secType, symbol=self.req.m_symbol, currency=self.req.m_currency, exchange=self.req.m_exchange, \
                   primaryExchange=self.req.m_primaryExch, endtime=self.req.endtime, duration=self.req.duration, \
                   bar_size=self.req.bar_size, what_to_show=self.req.what_to_show, use_rth=self.req.use_rth)

    def _error_handler(self, msg):
        print("error: ", msg)
        if not msg:
            print('disconnecting', self._conn.disconnect())



    def __init__(self, data_path='/docker_stocks/data', date_format='%Y%m%d %H:%M:%S', host='127.0.0.1', port=4001, client_id=None):
        self._data_path = data_path
        self._date_format = date_format
        self._next_valid_id = 1

        self._host = host
        self._port = port
        self._client_id = client_id






    def _historical_data_handler(self, msg):
        """
            Define historical data handler for IB - this will populate our pandas data frame
        """

        # print (msg.reqId, msg.date, msg.open, msg.close, msg.high, msg.low)
        if not 'finished' in str(msg.date):
            #print ("historical_data_handler: ", msg)
            try:
                self._s = ([datetime.strptime(msg.date, self._date_format),
                            msg.open, msg.high, msg.low, msg.close, msg.volume, 0])
            except:
                #for dates only with no time to str:
                self._s = ([datetime.strptime(msg.date, "%Y%m%d"),
                            msg.open, msg.high, msg.low, msg.close, msg.volume, 0])
            self._df.loc[len(self._df)] = self._s

        else:
            self._df.set_index('Date', inplace=True)



    def setAllWithKwArgs(self, **kwargs):
        #set all attributes to the kwargs to pass along
        for key, value in kwargs.items():
            setattr(self, key, value)

    def inner(self, sec_type, symbol, currency, exchange, primaryExchange, endtime, duration, bar_size, what_to_show, use_rth):
        print ("calling inner... setting up req.")
        self.req = Contract()
        self.req.m_secType = sec_type
        self.req.m_symbol = symbol
        self.req.m_currency = currency
        self.req.m_exchange = exchange
        self.primaryExch = primaryExchange
        self.req.endtime = endtime
        self.req.duration = duration
        self.req.bar_size = bar_size
        self.req.what_to_show = what_to_show
        self.req.use_rth = use_rth


        self._conn.reqHistoricalData(self._next_valid_id, self.req, endtime, duration,
                                     bar_size, what_to_show, use_rth, 1)


    def get_dataframe(self, sec_type, symbol, currency, exchange, primaryExchange, endtime, duration, bar_size, what_to_show, use_rth, timeoutsecs):



        # build filename
        self.filename = symbol + '_' + sec_type + '_' + exchange + '_' + currency + '_' + \
            endtime.replace(' ', '') + '_' + duration.replace(' ', '') + '_' + bar_size.replace(' ', '') + '_' + \
            what_to_show + '_' + str(use_rth) + '.csv'
        self.filename = self.filename.replace('/', '.')
        self.filename = self._data_path + '/' + self.filename
        print ("filename:  ", self.filename)


        try:

            # check if we have this cached
            self._df = pd.read_csv(self.filename,
                         parse_dates=True,
                         index_col=0)
            self._df.index = pd.to_datetime(self._df.index, format='%Y-%m-%d %H:%M:%S')
        except IOError:
            #set up connection:
            self._reset_data(self._host, self._port, self._client_id)

            # Not cached. Download it.
            # Establish a Contract object and the params for the request
            self.inner(sec_type, symbol, currency, exchange,
                  primaryExchange, endtime, duration, bar_size,
                  what_to_show, use_rth)

            # Make sure the connection doesn't get disconnected prior the response data return
            timeout = time() + timeoutsecs
            while self._conn.isConnected() and time() < timeout:
                #print(".", end="", flush=True)
                pass
            self._conn.disconnect()
            self._df.to_csv(self.filename)


        return self._df





if __name__ == '__main__':

    date_format = '%Y%m%d %H:%M:%S'

    downloader_kwargs = dict(
        data_path='../data',
        date_format=date_format,
        host='127.0.0.1',
        port=4001,
        client_id=992
    )
    downloader = IBDataCache(**downloader_kwargs)

    stock_kwargs = dict(
        sec_type='STK',
        symbol='AAPL',
        currency='USD',
        exchange='SMART',
        primaryExchange='NASDAQ',
        endtime=datetime(2018, 10, 26, 15, 59).strftime(date_format),
        duration='2 D',
        bar_size='30 mins',
        what_to_show='TRADES',
        use_rth=1
    )

    df = downloader.get_dataframe(**stock_kwargs)
    print(df)


    stock_kwargs = dict(
        sec_type='STK',
        symbol='MSFT',
        currency='USD',
        exchange='SMART',
        primaryExchange='NASDAQ',
        endtime=datetime(2018, 10, 26, 15, 59).strftime(date_format),
        duration='1 D',
        bar_size='1 min',
        what_to_show='TRADES',
        use_rth=1
    )

    df = downloader.get_dataframe(**stock_kwargs)

    print ("IBCacheData")

    print(df)

(As a side note, the int64 error mentioned above is due to the sleep statement from the original code terminating before the downloading is completed.)

Benoît -- have you tried runonce=False?

Hi,

Here's a wavelet bandpass filter indicator using the pywt (python wavelets) library. I haven't found the backtrader_contrib github so I'll just leave this here in case anyone wants to use it. GPL3. No support, warranties, either expressed or implied. No guarantee that this will have efficacy for your particular application.

With minimal (or no) changes this should support:

Haar (haar)
Daubechies (db)
Symlets (sym)
Coiflets (coif)
Biorthogonal (bior)
Reverse biorthogonal (rbio)
“Discrete” FIR approximation of Meyer wavelet (dmey)
Gaussian wavelets (gaus)
Mexican hat wavelet (mexh)
Morlet wavelet (morl)
Complex Gaussian wavelets (cgau)
Shannon wavelets (shan)
Frequency B-Spline wavelets (fbsp)
Complex Morlet wavelets (cmor)

Good luck,
B. Bradford

#!/usr/bin/env python
# -*- coding: utf-8; py-indent-offset:4 -*-
##############################################################################
#
# Copyright (C) 2017 Benjamin Bradford  (bigdavediode2 (at ) yahoo.com)
# Wavelet deconstruction to specified subband/level
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU 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 General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

import backtrader.indicators as btind
import backtrader as bt
import pywt

class bbDWT(bt.Indicator):

#   Approximation, detail coefficient lines if wanted:
    lines = ('cA2', 'cD2', 'cD1', )
    params = (('wavelettype', 'db2'),
             ('coefflevels', 5),
             ('coeffpick', 0), )




    def once(self, start, end):

        super(bbDWT, self).__init__()
        lstinputdata = self.data.array


        coeffs = pywt.wavedecn(lstinputdata, self.params.wavelettype, level=self.params.coefflevels)

#        Use if you wish to graph spectrum bands/coefficients:
#        coeffcA2, coeffcD2, coeffcD1 = coeffs
#        print(coeffs)
        print(len(coeffs))

        arr, coeff_slices = pywt.coeffs_to_array(coeffs)
        print(coeff_slices)

        # Zero out all but selected wavelet coeff level (avoids having to pad arrays):
        for i in range(0, len(coeff_slices)):
            if i != self.params.coeffpick:
                #Remove the "slice" portion of the returned slices up to the bracket to run as a slice command:
                s1 = str(coeff_slices[i])
                s1 = s1[s1.find('slice'):s1.find(')')+1]
                execstring = 's1 = ' + s1
                print(execstring)
                exec(execstring)
                arr[s1] = 0

        print(arr)


        coeffs_from_array = pywt.array_to_coeffs(arr, coeff_slices)

        dummy_array = pywt.waverecn(coeffs_from_array, self.params.wavelettype)

        self.lines.cA2.array = dummy_array

@incabot Or perhaps convert your Pandas dataframe datetime column to an actual datetime index rather than a string.

For anyone who needs it, here is an indicator that calculates and charts implied volatility. I used GBS/Bjerksund-Stensland 2002 routines from the incredibly helpful Davis Edwards, author of the book "Risk Management in Trading."

It calculates implied volatility for American options only but could easily be adjusted to support European options. It is currently set to download the underlying price (data0) and the options data (data1) from IB.

In lieu of IB's IV data this allows you to adjust the underlying cost of carry/dividend rate/risk free rate using overnight index swaps or whatever you prefer. Risk free rate for American options should be taken from Overnight Indexed Swaps -- e.g. the Federal Funds rate which is published daily by the Federal Reserve in the US. Overnight rates include EONIA (EUR), SONIA (GBP), CHOIS (CHF), and TONAR (JPY). Recommended to not use LIBOR nor the T-Bill rate.

If it's not working for you, check that your option string uses the IB format: e.g. 'AAPL-20180119-SMART-USD-170-CALL' and that you are not trying to download a weekly that cannot be routed through SMART. Also check your expiry dates at the IB website, as they might be slightly off from what might be expected.

Here is the indicator:

'''
Author:  Davis Edwards ported to backtrader by Benjamin Bradford 

MIT License

Copyright (c) 2017

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
'''


import math
import numpy as np
from scipy.stats import norm
from scipy.stats import mvn
import datetime
import backtrader as bt
import backtrader.indicators as btind


# This class contains the limits on inputs for GBS models
# It is not intended to be part of this module's public interface
class _GBS_Limits:
    # An GBS model will return an error if an out-of-bound input is input
    MAX32 = 2147483248.0

    MIN_T = 1.0 / 1000.0  # requires some time left before expiration
    MIN_X = 0.01
    MIN_FS = 0.01

    # Volatility smaller than 0.5% causes American Options calculations
    # to fail (Number to large errors).
    # GBS() should be OK with any positive number. Since vols less
    # than 0.5% are expected to be extremely rare, and most likely bad inputs,
    # _gbs() is assigned this limit too
    MIN_V = 0.005

    MAX_T = 100
    MAX_X = MAX32
    MAX_FS = MAX32

    # Asian Option limits
    # maximum TA is time to expiration for the option
    MIN_TA = 0

    # This model will work with higher values for b, r, and V. However, such values are extremely uncommon.
    # To catch some common errors, interest rates and volatility is capped to 100%
    # This reason for 1 (100%) is mostly to cause the library to throw an exceptions
    # if a value like 15% is entered as 15 rather than 0.15)
    MIN_b = -1
    MIN_r = -1

    MAX_b = 1
    MAX_r = 1
    MAX_V = 1


# The primary class for calculating Generalized Black Scholes option prices and deltas
# It is not intended to be part of this module's public interface

# Inputs: option_type = "p" or "c", fs = price of underlying, x = strike, t = time to expiration, r = risk free rate
#         b = cost of carry, v = implied volatility
# Outputs: value, delta, gamma, theta, vega, rho
def _gbs(option_type, fs, x, t, r, b, v):
#    _debug("Debugging Information: _gbs()")
    # -----------
    # Test Inputs (throwing an exception on failure)
    _gbs_test_inputs(option_type, fs, x, t, r, b, v)

    # -----------
    # Create preliminary calculations
    t__sqrt = math.sqrt(t)
    d1 = (math.log(fs / x) + (b + (v * v) / 2) * t) / (v * t__sqrt)
    d2 = d1 - v * t__sqrt

    if option_type == "c":
        # it's a call
#        _debug("     Call Option")
        value = fs * math.exp((b - r) * t) * norm.cdf(d1) - x * math.exp(-r * t) * norm.cdf(d2)
        delta = math.exp((b - r) * t) * norm.cdf(d1)
        gamma = math.exp((b - r) * t) * norm.pdf(d1) / (fs * v * t__sqrt)
        theta = -(fs * v * math.exp((b - r) * t) * norm.pdf(d1)) / (2 * t__sqrt) - (b - r) * fs * math.exp(
            (b - r) * t) * norm.cdf(d1) - r * x * math.exp(-r * t) * norm.cdf(d2)
        vega = math.exp((b - r) * t) * fs * t__sqrt * norm.pdf(d1)
        rho = x * t * math.exp(-r * t) * norm.cdf(d2)
    else:
        # it's a put
#        _debug("     Put Option")
        value = x * math.exp(-r * t) * norm.cdf(-d2) - (fs * math.exp((b - r) * t) * norm.cdf(-d1))
        delta = -math.exp((b - r) * t) * norm.cdf(-d1)
        gamma = math.exp((b - r) * t) * norm.pdf(d1) / (fs * v * t__sqrt)
        theta = -(fs * v * math.exp((b - r) * t) * norm.pdf(d1)) / (2 * t__sqrt) + (b - r) * fs * math.exp(
            (b - r) * t) * norm.cdf(-d1) + r * x * math.exp(-r * t) * norm.cdf(-d2)
        vega = math.exp((b - r) * t) * fs * t__sqrt * norm.pdf(d1)
        rho = -x * t * math.exp(-r * t) * norm.cdf(-d2)

#    _debug("     d1= {0}\n     d2 = {1}".format(d1, d2))
#    _debug("     delta = {0}\n     gamma = {1}\n     theta = {2}\n     vega = {3}\n     rho={4}".format(delta, gamma,
#                                                                                                        theta, vega,
#                                                                                                        rho))

    return value, delta, gamma, theta, vega, rho

# ------------------------------
# This function verifies that the Call/Put indicator is correctly entered
def _test_option_type(option_type):
    if (option_type != "c") and (option_type != "p"):
        raise GBS_InputError("Invalid Input option_type ({0}). Acceptable value are: c, p".format(option_type))

# ------------------------------
# This function makes sure inputs are OK
# It throws an exception if there is a failure
def _gbs_test_inputs(option_type, fs, x, t, r, b, v):
    # -----------
    # Test inputs are reasonable
    _test_option_type(option_type)

    if (x < _GBS_Limits.MIN_X) or (x > _GBS_Limits.MAX_X):
        raise GBS_InputError(
            "Invalid Input Strike Price (X). Acceptable range for inputs is {1} to {2}".format(x, _GBS_Limits.MIN_X,
                                                                                               _GBS_Limits.MAX_X))

    if (fs < _GBS_Limits.MIN_FS) or (fs > _GBS_Limits.MAX_FS):
        raise GBS_InputError(
            "Invalid Input Forward/Spot Price (FS). Acceptable range for inputs is {1} to {2}".format(fs,
                                                                                                      _GBS_Limits.MIN_FS,
                                                                                                      _GBS_Limits.MAX_FS))

    if (t < _GBS_Limits.MIN_T) or (t > _GBS_Limits.MAX_T):
        raise GBS_InputError(
            "Invalid Input Time (T = {0}). Acceptable range for inputs is {1} to {2}".format(t, _GBS_Limits.MIN_T,
                                                                                             _GBS_Limits.MAX_T))

    if (b < _GBS_Limits.MIN_b) or (b > _GBS_Limits.MAX_b):
        raise GBS_InputError(
            "Invalid Input Cost of Carry (b = {0}). Acceptable range for inputs is {1} to {2}".format(b,
                                                                                                      _GBS_Limits.MIN_b,
                                                                                                      _GBS_Limits.MAX_b))

    if (r < _GBS_Limits.MIN_r) or (r > _GBS_Limits.MAX_r):
        raise GBS_InputError(
            "Invalid Input Risk Free Rate (r = {0}). Acceptable range for inputs is {1} to {2}".format(r,
                                                                                                       _GBS_Limits.MIN_r,
                                                                                                       _GBS_Limits.MAX_r))

    if (v < _GBS_Limits.MIN_V) or (v > _GBS_Limits.MAX_V):
        raise GBS_InputError(
            "Invalid Input Implied Volatility (V = {0}). Acceptable range for inputs is {1} to {2}".format(v,
                                                                                                           _GBS_Limits.MIN_V,
                                                                                                           _GBS_Limits.MAX_V))
# This class defines the Exception that gets thrown when invalid input is placed into the GBS function
class GBS_InputError(Exception):
    def __init__(self, mismatch):
        Exception.__init__(self, mismatch)

# This class defines the Exception that gets thrown when there is a calculation error
class GBS_CalculationError(Exception):
    def __init__(self, mismatch):
        Exception.__init__(self, mismatch)

class bbImpVol(bt.Indicator):

    lines = ('impvolplot', )

    # Riskfreerate for American options should be taken from
    # Overnight Indexed Swaps -- e.g. the Federal Funds rate which
    # is published daily by the Federal Reserve in the US. Overnight
    # rates include EONIA (EUR), SONIA (GBP), CHOIS (CHF), and TONAR (JPY).
    # Do not use LIBOR nor the T-Bill rate.
    params = (('optstring', 'AAPL-20180119-SMART-USD-170-CALL'),
              ('riskfreerate', 0.0142),
              ('dividends', 0.0143))    # Dividend yield p.a. for AAPL


    # ----------
    # Implied Volatility Calculations:
    # Inputs (not all functions use all inputs)
    #      fs = forward/spot price
    #      x = Strike
    #      t = Time (in years)
    #      r = risk free rate
    #      b = cost of carry
    #      cp = Call or Put price
    #      precision = (optional) precision at stopping point
    #      max_steps = (optional) maximum number of steps

    # ----------
    # Approximate Implied Volatility
    #
    # This function is used to choose a starting point for the
    # search functions (Newton and bisection searches).
    # Brenner & Subrahmanyam (1988), Feinstein (1988)


    # -----------
    # Generalized American Option Pricer
    # This is a wrapper to check inputs and route to the current "best" American option model
    def _american_option(self, option_type, fs, x, t, r, b, v):
        # -----------
        # Test Inputs (throwing an exception on failure)
#        _debug("Debugging Information: _american_option()")
        _gbs_test_inputs(option_type, fs, x, t, r, b, v)

        # -----------
        if option_type == "c":
            # Call Option
#            _debug("     Call Option")
            return self._bjerksund_stensland_2002(fs, x, t, r, b, v)
        else:
            # Put Option
#            _debug("     Put Option")

            # Using the put-call transformation: P(X, FS, T, r, b, V) = C(FS, X, T, -b, r-b, V)
            # WARNING - When reconciling this code back to the B&S paper, the order of variables is different

            put__x = fs
            put_fs = x
            put_b = -b
            put_r = r - b

            # pass updated values into the Call Valuation formula
            return self._bjerksund_stensland_2002(put_fs, put__x, t, put_r, put_b, v)

    # -----------
    # American Call Option (Bjerksund Stensland 2002 approximation)
    def _bjerksund_stensland_2002(self, fs, x, t, r, b, v):
        # -----------
        # initialize output
        # using GBS greeks (TO DO: update greek calculations)
        my_output = _gbs("c", fs, x, t, r, b, v)

        e_value = my_output[0]
        delta = my_output[1]
        gamma = my_output[2]
        theta = my_output[3]
        vega = my_output[4]
        rho = my_output[5]

        # debugging for calculations
#        _debug("-----")
#        _debug("Debug Information: _Bjerksund_Stensland_2002())")

        # If b >= r, it is never optimal to exercise before maturity
        # so we can return the GBS value
        if b >= r:
            #print ("     Returning GBS value")
            #print e_value, delta, gamma, theta, vega, rho
            return e_value, delta, gamma, theta, vega, rho

        # -----------
        # Create preliminary calculations
        print ("Returning Bjerk/Sten value")
        v2 = v ** 2
        t1 = 0.5 * (math.sqrt(5) - 1) * t
        t2 = t

        beta_inside = ((b / v2 - 0.5) ** 2) + 2 * r / v2
        # forcing the inside of the sqrt to be a positive number
        beta_inside = abs(beta_inside)
        beta = (0.5 - b / v2) + math.sqrt(beta_inside)
        b_infinity = (beta / (beta - 1)) * x
        b_zero = max(x, (r / (r - b)) * x)

        h1 = -(b * t1 + 2 * v * math.sqrt(t1)) * ((x ** 2) / ((b_infinity - b_zero) * b_zero))
        h2 = -(b * t2 + 2 * v * math.sqrt(t2)) * ((x ** 2) / ((b_infinity - b_zero) * b_zero))

        i1 = b_zero + (b_infinity - b_zero) * (1 - math.exp(h1))
        i2 = b_zero + (b_infinity - b_zero) * (1 - math.exp(h2))

        alpha1 = (i1 - x) * (i1 ** (-beta))
        alpha2 = (i2 - x) * (i2 ** (-beta))

        # debugging for calculations
#        _debug("     t1 = {0}".format(t1))
#        _debug("     beta = {0}".format(beta))
#        _debug("     b_infinity = {0}".format(b_infinity))
#        _debug("     b_zero = {0}".format(b_zero))
#        _debug("     h1 = {0}".format(h1))
#        _debug("     h2 = {0}".format(h2))
#        _debug("     i1 = {0}".format(i1))
#        _debug("     i2 = {0}".format(i2))
#        _debug("     alpha1 = {0}".format(alpha1))
#        _debug("     alpha2 = {0}".format(alpha2))

        # check for immediate exercise
        if fs >= i2:
            value = fs - x
        else:
            # Perform the main calculation
            value = (alpha2 * (fs ** beta)
                     - alpha2 * self._phi(fs, t1, beta, i2, i2, r, b, v)
                     + self._phi(fs, t1, 1, i2, i2, r, b, v)
                     - self._phi(fs, t1, 1, i1, i2, r, b, v)
                     - x * self._phi(fs, t1, 0, i2, i2, r, b, v)
                     + x * self._phi(fs, t1, 0, i1, i2, r, b, v)
                     + alpha1 * self._phi(fs, t1, beta, i1, i2, r, b, v)
                     - alpha1 * self._psi(fs, t2, beta, i1, i2, i1, t1, r, b, v)
                     + self._psi(fs, t2, 1, i1, i2, i1, t1, r, b, v)
                     - self._psi(fs, t2, 1, x, i2, i1, t1, r, b, v)
                     - x * self._psi(fs, t2, 0, i1, i2, i1, t1, r, b, v)
                     + x * self._psi(fs, t2, 0, x, i2, i1, t1, r, b, v))

        # in boundary conditions, this approximation can break down
        # Make sure option value is greater than or equal to European value
        value = max(value, e_value)

        # -----------
        # Return Data
        return value, delta, gamma, theta, vega, rho


    # ---------------------------
    # American Option Intermediate Calculations

    # -----------
    # The Psi() function used by _Bjerksund_Stensland_2002 model
    def _psi(self, fs, t2, gamma, h, i2, i1, t1, r, b, v):
        vsqrt_t1 = v * math.sqrt(t1)
        vsqrt_t2 = v * math.sqrt(t2)

        bgamma_t1 = (b + (gamma - 0.5) * (v ** 2)) * t1
        bgamma_t2 = (b + (gamma - 0.5) * (v ** 2)) * t2

        d1 = (math.log(fs / i1) + bgamma_t1) / vsqrt_t1
        d3 = (math.log(fs / i1) - bgamma_t1) / vsqrt_t1

        d2 = (math.log((i2 ** 2) / (fs * i1)) + bgamma_t1) / vsqrt_t1
        d4 = (math.log((i2 ** 2) / (fs * i1)) - bgamma_t1) / vsqrt_t1

        e1 = (math.log(fs / h) + bgamma_t2) / vsqrt_t2
        e2 = (math.log((i2 ** 2) / (fs * h)) + bgamma_t2) / vsqrt_t2
        e3 = (math.log((i1 ** 2) / (fs * h)) + bgamma_t2) / vsqrt_t2
        e4 = (math.log((fs * (i1 ** 2)) / (h * (i2 ** 2))) + bgamma_t2) / vsqrt_t2

        tau = math.sqrt(t1 / t2)
        lambda1 = (-r + gamma * b + 0.5 * gamma * (gamma - 1) * (v ** 2))
        kappa = (2 * b) / (v ** 2) + (2 * gamma - 1)

        psi = math.exp(lambda1 * t2) * (fs ** gamma) * (self._cbnd(-d1, -e1, tau)
                                                        - ((i2 / fs) ** kappa) * self._cbnd(-d2, -e2, tau)
                                                        - ((i1 / fs) ** kappa) * self._cbnd(-d3, -e3, -tau)
                                                        + ((i1 / i2) ** kappa) * self._cbnd(-d4, -e4, -tau))
        return psi


    # -----------
    # The Phi() function used by _Bjerksund_Stensland_2002 model and the _Bjerksund_Stensland_1993 model
    def _phi(self, fs, t, gamma, h, i, r, b, v):
        d1 = -(math.log(fs / h) + (b + (gamma - 0.5) * (v ** 2)) * t) / (v * math.sqrt(t))
        d2 = d1 - 2 * math.log(i / fs) / (v * math.sqrt(t))

        lambda1 = (-r + gamma * b + 0.5 * gamma * (gamma - 1) * (v ** 2))
        kappa = (2 * b) / (v ** 2) + (2 * gamma - 1)

        phi = math.exp(lambda1 * t) * (fs ** gamma) * (norm.cdf(d1) - ((i / fs) ** kappa) * norm.cdf(d2))

#        _debug("-----")
#        _debug("Debug info for: _phi()")
#        _debug("    d1={0}".format(d1))
#        _debug("    d2={0}".format(d2))
#        _debug("    lambda={0}".format(lambda1))
#        _debug("    kappa={0}".format(kappa))
#        _debug("    phi={0}".format(phi))
        return phi


    # -----------
    # Cumulative Bivariate Normal Distribution
    # Primarily called by Psi() function, part of the _Bjerksund_Stensland_2002 model
    def _cbnd(self, a, b, rho):
        # This distribution uses the Genz multi-variate normal distribution
        # code found as part of the standard SciPy distribution
        lower = np.array([0, 0])
        upper = np.array([a, b])
        infin = np.array([0, 0])
        correl = rho
        error, value, inform = mvn.mvndst(lower, upper, infin, correl)
        return value


    def _approx_implied_vol(self, option_type, fs, x, t, r, b, cp):
        self._test_option_type(option_type)

        ebrt = math.exp((b - r) * t)
        ert = math.exp(-r * t)

        a = math.sqrt(2 * math.pi) / (fs * ebrt + x * ert)

        if option_type == "c":
            payoff = fs * ebrt - x * ert
        else:
            payoff = x * ert - fs * ebrt

        b = cp - payoff / 2
        c = (payoff ** 2) / math.pi

        v = (a * (b + math.sqrt(b ** 2 + c))) / math.sqrt(t)

        return v

    # ------------------------------
    # This function verifies that the Call/Put indicator is correctly entered
    def _test_option_type(self, option_type):
        if (option_type != "c") and (option_type != "p"):
            raise GBS_InputError("Invalid Input option_type ({0}). Acceptable value are: c, p".format(option_type))


    # ------------------------------
    # This function makes sure inputs are OK
    # It throws an exception if there is a failure
    def _gbs_test_inputs(option_type, fs, x, t, r, b, v):
        # -----------
        # Test inputs are reasonable
        _test_option_type(option_type)

        if (x < _GBS_Limits.MIN_X) or (x > _GBS_Limits.MAX_X):
            raise GBS_InputError(
                "Invalid Input Strike Price (X). Acceptable range for inputs is {1} to {2}".format(x, _GBS_Limits.MIN_X,
                                                                                                   _GBS_Limits.MAX_X))

        if (fs < _GBS_Limits.MIN_FS) or (fs > _GBS_Limits.MAX_FS):
            raise GBS_InputError(
                "Invalid Input Forward/Spot Price (FS). Acceptable range for inputs is {1} to {2}".format(fs,
                                                                                                          _GBS_Limits.MIN_FS,
                                                                                                          _GBS_Limits.MAX_FS))

        if (t < _GBS_Limits.MIN_T) or (t > _GBS_Limits.MAX_T):
            raise GBS_InputError(
                "Invalid Input Time (T = {0}). Acceptable range for inputs is {1} to {2}".format(t, _GBS_Limits.MIN_T,
                                                                                                 _GBS_Limits.MAX_T))

        if (b < _GBS_Limits.MIN_b) or (b > _GBS_Limits.MAX_b):
            raise GBS_InputError(
                "Invalid Input Cost of Carry (b = {0}). Acceptable range for inputs is {1} to {2}".format(b,
                                                                                                          _GBS_Limits.MIN_b,
                                                                                                          _GBS_Limits.MAX_b))

        if (r < _GBS_Limits.MIN_r) or (r > _GBS_Limits.MAX_r):
            raise GBS_InputError(
                "Invalid Input Risk Free Rate (r = {0}). Acceptable range for inputs is {1} to {2}".format(r,
                                                                                                           _GBS_Limits.MIN_r,
                                                                                                           _GBS_Limits.MAX_r))

        if (v < _GBS_Limits.MIN_V) or (v > _GBS_Limits.MAX_V):
            raise GBS_InputError(
                "Invalid Input Implied Volatility (V = {0}). Acceptable range for inputs is {1} to {2}".format(v,
                                                                                                               _GBS_Limits.MIN_V,
                                                                                                               _GBS_Limits.MAX_V))

    # ----------
    # Find the Implied Volatility of an European (GBS) Option given a price
    # using Newton-Raphson method for greater speed since Vega is available

    def _gbs_implied_vol(option_type, fs, x, t, r, b, cp, precision=.00001, max_steps=100):
        return _newton_implied_vol(_gbs, option_type, x, fs, t, b, r, cp, precision, max_steps)
    # ----------
    # Find the Implied Volatility of an American Option given a price
    # Using bisection method since Vega is difficult to estimate for Americans
    def _american_implied_vol(self, option_type, fs, x, t, r, b, cp, precision=.00001, max_steps=100):
        return self._bisection_implied_vol(self._american_option, option_type, fs, x, t, r, b, cp, precision, max_steps)


    # ----------
    # Calculate Implied Volatility with a Newton Raphson search
    def _newton_implied_vol(val_fn, option_type, x, fs, t, b, r, cp, precision=.00001, max_steps=100):
        # make sure a valid option type was entered
        _test_option_type(option_type)

        # Estimate starting Vol, making sure it is allowable range
        v = _approx_implied_vol(option_type, fs, x, t, r, b, cp)
        v = max(_GBS_Limits.MIN_V, v)
        v = min(_GBS_Limits.MAX_V, v)

        # Calculate the value at the estimated vol
        value, delta, gamma, theta, vega, rho = val_fn(option_type, fs, x, t, r, b, v)
        min_diff = abs(cp - value)

#        _debug("-----")
#        _debug("Debug info for: _Newton_ImpliedVol()")
#        _debug("    Vinitial={0}".format(v))

        # Newton-Raphson Search
        countr = 0
        while precision <= abs(cp - value) <= min_diff and countr < max_steps:

            v = v - (value - cp) / vega
            if (v > _GBS_Limits.MAX_V) or (v < _GBS_Limits.MIN_V):
#                _debug("    Volatility out of bounds")
                break

            value, delta, gamma, theta, vega, rho = val_fn(option_type, fs, x, t, r, b, v)
            min_diff = min(abs(cp - value), min_diff)

            # keep track of how many loops
            countr += 1
#            _debug("     IVOL STEP {0}. v={1}".format(countr, v))

        # check if function converged and return a value
        if abs(cp - value) < precision:
            # the search function converged
            return v
        else:
            # if the search function didn't converge, try a bisection search
            return _bisection_implied_vol(val_fn, option_type, fs, x, t, r, b, cp, precision, max_steps)


    # ----------
    # Find the Implied Volatility using a Bisection search (American options)
    def _bisection_implied_vol(self, val_fn, option_type, fs, x, t, r, b, cp, precision=.00001, max_steps=100):


        # Estimate Upper and Lower bounds on volatility
        # Assume American Implied vol is within +/- 50% of the GBS Implied Vol
        v_mid = self._approx_implied_vol(option_type, fs, x, t, r, b, cp)

        if (v_mid <= _GBS_Limits.MIN_V) or (v_mid >= _GBS_Limits.MAX_V):
            # if the volatility estimate is out of bounds, search entire allowed vol space
            v_low = _GBS_Limits.MIN_V
            v_high = _GBS_Limits.MAX_V
            v_mid = (v_low + v_high) / 2
        else:
            # reduce the size of the vol space
            v_low = max(_GBS_Limits.MIN_V, v_mid * .5)
            v_high = min(_GBS_Limits.MAX_V, v_mid * 1.5)

        # Estimate the high/low bounds on price
        cp_mid = val_fn(option_type, fs, x, t, r, b, v_mid)[0]

        # initialize bisection loop
        current_step = 0
        diff = abs(cp - cp_mid)

#        _debug("     American IVOL starting conditions: CP={0} cp_mid={1}".format(cp, cp_mid))
#        _debug("     IVOL {0}. V[{1},{2},{3}]".format(current_step, v_low, v_mid, v_high))

        # Keep bisection volatility until correct price is found
        while (diff > precision) and (current_step < max_steps):
            current_step += 1

            # Cut the search area in half
            if cp_mid < cp:
                v_low = v_mid
            else:
                v_high = v_mid

            cp_low = val_fn(option_type, fs, x, t, r, b, v_low)[0]
            cp_high = val_fn(option_type, fs, x, t, r, b, v_high)[0]

            v_mid = v_low + (cp - cp_low) * (v_high - v_low) / (cp_high - cp_low)
            v_mid = max(_GBS_Limits.MIN_V, v_mid)  # enforce high/low bounds
            v_mid = min(_GBS_Limits.MAX_V, v_mid)  # enforce high/low bounds

            cp_mid = val_fn(option_type, fs, x, t, r, b, v_mid)[0]
            diff = abs(cp - cp_mid)

#            _debug("     IVOL {0}. V[{1},{2},{3}]".format(current_step, v_low, v_mid, v_high))

        # return output
        if abs(cp - cp_mid) < precision:
            return v_mid
        else:
            raise GBS_CalculationError(
                "Implied Vol did not converge. Best Guess={0}, Price diff={1}, Required Precision={2}".format(v_mid,
                                                                                                              diff,
                                                                                                              precision))
# Public Interface for implied Volatility Functions
    # Inputs:
    #    option_type = "p" or "c"
    #    fs          = price of underlying
    #    x           = strike
    #    t           = time to expiration
    #    v           = implied volatility
    #    r           = risk free rate
    #    q           = dividend payment
    #    b           = cost of carry
    # Outputs:
    #    value       = price of the option
    #    delta       = first derivative of value with respect to price of underlying
    #    gamma       = second derivative of value w.r.t price of underlying
    #    theta       = first derivative of value w.r.t. time to expiration
    #    vega        = first derivative of value w.r.t. implied volatility
    #    rho         = first derivative of value w.r.t. risk free rates

    def euro_implied_vol(option_type, fs, x, t, r, q, cp):
        b = r - q
        return _gbs_implied_vol(option_type, fs, x, t, r, b, cp)


    def euro_implied_vol_76(option_type, fs, x, t, r, cp):
    # for European commodity options
        b = 0
        return _gbs_implied_vol(option_type, fs, x, t, r, b, cp)


    def amer_implied_vol(self, option_type, fs, x, t, r, q, cp):
        b = r - q
        return self._american_implied_vol(option_type, fs, x, t, r, b, cp)


    def amer_implied_vol_76(option_type, fs, x, t, r, cp):
    # for American commodity options
        b = 0
        return _american_implied_vol(option_type, fs, x, t, r, b, cp)


    def __init__(self):
        #Set the option information, strike, expiry, etc.
        cpinf = self.params.optstring.split('-')
        #          print cpinf
        #          print cpinf[1]

        self.expdate = datetime.datetime.strptime(cpinf[1], "%Y%m%d").date()
        # opt type, underlying price, strike, time to expiry, riskfreerate, cost of carry, call or put price:
        self.sopttype = cpinf[5][:1].lower()
        self.doptstrike = float(cpinf[4])
        print self.expdate, self.sopttype, self.doptstrike

    def next(self):
        print (self.params.optstring)

        # print type(self.params.optiontype), type(self.params.underlyingprice), type(self.params.strike), type(self.params.timeexp), type(self.params.riskfreerate), type(self.params.costcarry), type(self.params.cpprice)

        self.lines.impvolplot.subplot=True
        print "--------"

        #         time to expiration is in years, and must be a float:
        self.time_d = self.expdate - self.data.datetime.date(ago=0)
        print self.time_d, self.expdate, self.data.datetime.date(ago=0)
        timeexp = float(self.time_d.days) / 365
        print timeexp

        #        print "try except: ", self.datas[1].close[-1]
        #        try:
        self.callputprice = float(self.datas[1].close[0])
        #        except:
        #            self.callputprice = float(self.datas[1].close[-1])

        self.underlyingprice = float(self.datas[0].close[0])

        kwargs = dict(
            option_type=self.sopttype,
            fs=self.underlyingprice,
            x=self.doptstrike,
            t=timeexp,
            r=self.params.riskfreerate,
            q=self.params.dividends,
            cp=self.callputprice, )
        print "imp vol is: ", bbImpVol().amer_implied_vol(**kwargs)
        self.lines.impvolplot[0] = bbImpVol().amer_implied_vol(**kwargs)
        print "self lines imp vol is: ", self.lines.impvolplot[0]


        print "End of bbImpliedVolatility"
        print "----------------"

Fixed a real time bar delay bug and tightened the code:

    def __init__(self):

        cpy = copy.copy(self.data)
        tmp = bt.If(self.data(0) == self.data(-1), cpy(-1) + 0.000001, self.data(0))
        self.lines.zigzag_peak = bt.If(bt.And(tmp(-1)>tmp(-2), tmp(-1)>tmp(0)), self.data(-1), float('nan'))
        cpy = copy.copy(self.data)
        tmp = bt.If(self.data(0) == self.data(-1), cpy(-1) - 0.000001, self.data(0))
        self.lines.zigzag_valley = bt.If(bt.And(tmp(-1) < tmp(-2), tmp(-1) < tmp(0)), self.data(-1), float('nan'))
        self.lines.zigzag = bt.If(self.zigzag_peak, self.zigzag_peak, bt.If(self.zigzag_valley, self.zigzag_valley, float('nan')))

Very nice. Someday I'd like to meet Daniel R. in person and buy him a beverage.

Vlad, I'm limited in time and scope but would like to try to contribute as an admin where I can: https://github.com/bigdavediode

Hi Kevin --

Another option for you is to use the IB data cache code that brettelliot wrote at https://community.backtrader.com/topic/1007/how-to-cached-data-from-ib-and-use-it-later/2

It exposes the IB connection so that the ib api parameter for RTH is available (use_rth).

Also gets you a nice cache and a bit of a speedup as well.

B.

@incabot Or perhaps convert your Pandas dataframe datetime column to an actual datetime index rather than a string.

@incabot Perhaps on the SQLLite3 side: PRAGMA table_info(stock_price);

And see if the columns are strings

Here's the code I previously posted to do exactly that:

https://community.backtrader.com/topic/785/calculate-an-option-s-implied-volatility-using-gbs-bjerksund-stensland-2002

@Armaan-Syed

I don't have example code to do that in this particular case because I don't reset the VWAP each day, I just use a running tally.

You may wish to consult this thread which will show you that getting datetimes in an indicator can be a challenge: https://community.backtrader.com/topic/489/getting-the-current-date-in-a-strategy-s-next-method/2

I have done it (iirc) by passing the datafeed separately into the indicator as datas[0] and the line that I'm operating the indicator on as datas[1].

So in my indicator I can filter (or in your case you can reset the VWAP) within the indicator based on the date:

curdt = self.datas[0].datetime.date(ago=0)
curtm = self.datas[0].datetime.time(ago=0)

And my line calculation input is in datas[1]. And I call passing both the data line and the input line that I'm working on:

x = DistPercentileSampled(d, self.inds[d]['highdiffs'], period=780, raw=True, bar_size=30, matchwkdy=False, inputdesc='SELL_disthighdiffs', subplot = False)

d is datas[0] and self.inds[d]['highdiffs'] is datas[1] within the indicator.

@Armaan-Syed This routine doesn't reset at the start of a trading day and continuously adjusts based on whatever period you supply. If you wish to have it reset at the start of each day you would need to access the datetime in the indicator, which might be problematic. An easier solution might just be to change typprice and cumtypprice to lines and just calculate a temporary daily running volume.

Or I suppose you could pass the data feed itself along with the actual data line to the indicator so you can recover the datetime information from the feed within your indicator to reset your running daily volume calculation at the start of day. I've passed those that way in the past with several indicators.

Having said that, I suggest that a running daily VWAP is more useful than one that resets daily. But I'm open to arguments.

@Martin-Bouhier said in ZigZag based on Close Prices:

Input("Reversal",0,100,5);

Hi Martin, I would suspect that you're either not using the same input data (ie. a moving average or smoother that is smaller and finer in the Metastock formula) or that the "percentage" input involves specifying a retracement perhaps before it counts as a peak or valley in that zig zag routine. Also I'm curious to know what the 5 in the "reversal" params specifies.

I don't like that python version. If there are repeated values at a peak, for example such as 200, 201, 201, 199, then it's possible for the code to entirely miss a peak.

Good luck.

Fixed a real time bar delay bug and tightened the code:

    def __init__(self):

        cpy = copy.copy(self.data)
        tmp = bt.If(self.data(0) == self.data(-1), cpy(-1) + 0.000001, self.data(0))
        self.lines.zigzag_peak = bt.If(bt.And(tmp(-1)>tmp(-2), tmp(-1)>tmp(0)), self.data(-1), float('nan'))
        cpy = copy.copy(self.data)
        tmp = bt.If(self.data(0) == self.data(-1), cpy(-1) - 0.000001, self.data(0))
        self.lines.zigzag_valley = bt.If(bt.And(tmp(-1) < tmp(-2), tmp(-1) < tmp(0)), self.data(-1), float('nan'))
        self.lines.zigzag = bt.If(self.zigzag_peak, self.zigzag_peak, bt.If(self.zigzag_valley, self.zigzag_valley, float('nan')))

Just to mention if you're attempting to use this live it obviously looks ahead one bar into the future which isn't allowed. If you're using it for preloaded historical data in a batch it works for identifying the exact turning point rather than one bar late. But for live data you'll need to move it one bar back by subtracting one from every bar reference as below. I added a couple of del lines just to prompt a bit of memory cleanup but they are entirely optional.

And yeah, I should use camelcase for the class name. Added a couple of comments to help with anyone who wants to understand the code.

   def __init__(self):

        #Make a copy
        tmp = copy.copy(self.data)
        #Peak shift
        tmp = bt.If(self.data(-1) == self.data(-2), tmp(-2) + 0.000001, self.data(-1))
        #Find peaks
        self.zigzag_peak = bt.If(bt.And(tmp(-1)>tmp(-2), tmp(-1)>tmp(0)), self.data(-1), float('nan'))
        del tmp
        tmp = copy.copy(self.data)
        #valley shift
        tmp = bt.If(self.data(-1) == self.data(-2), tmp(-2) - 0.000001, self.data(-1))
        #Find valleys:
        self.zigzag_valley = bt.If(bt.And(tmp(-1) < tmp(-2), tmp(-1) < tmp(0)), self.data(-1), float('nan'))
        self.lines.zigzag = bt.If(self.zigzag_peak, self.zigzag_peak,
                                  bt.If(self.zigzag_valley, self.zigzag_valley, float('nan')))
        self.lines.zigzag_peak = self.zigzag_peak * (1 + self.p.plotdistance)
        self.lines.zigzag_valley = self.zigzag_valley * (1 - self.p.plotdistance)
        del tmp

I realize this chart isn't the best illustration as the SMA wiggles and results in some doubling up of peak/valley indications which makes it difficult to see if it's showing a high or low at points.

@amulya-agarwal

Use compression, as above, compression 1 for 1 minute bars if you wish, as below:

        timeframe=bt.TimeFrame.Minutes,
        compression=barsize,
        useRTH=False,
        #outsideRth=True,    fill orders outside rth -- use for orders, not here...
        sessionstart=datetime.time(9, 30, 0),
        sessionend=datetime.time(16, 30, 0),

...and then your adddata or resample line using these arguments.