Python scipy 模块，exp() 实例源码

def devi(yy, eta):
    deveta = yy*eta - scipy.exp(eta)
    devy = yy*scipy.log(yy) - yy
    devy[yy == 0] = 0
    result = 2*(devy - deveta)
    return(result)

def rbf_kernel_pca(X, gamma, n_components):
    sq_dists = pdist(X, 'sqeuclidean')
    mat_sq_dists = squareform(sq_dists)
    K = exp(-gamma * mat_sq_dists)

    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)

    eigenvalues, eigenvectors = eigh(K)

    alphas = np.column_stack((
        eigenvectors[:, -i] for i in range(1, n_components+1)
    ))
    lambdas = [eigenvalues[-i] for i in range(1, n_components+1)]

    return alphas, lambdas

def KMV_f(E,D,T,r,sigmaE):
    n=10000
    m=2000
    diffOld=1e6     # a very big number
    for i in sp.arange(1,10):
        for j in sp.arange(1,m):
            A=E+D/2+i*D/n
            sigmaA=0.05+j*(1.0-0.001)/m
            d1 = (log(A/D)+(r+sigmaA*sigmaA/2.)*T)/(sigmaA*sqrt(T))
            d2 = d1-sigmaA*sqrt(T)
            diff4E= (A*N(d1)-D*exp(-r*T)*N(d2)-E)/A  # scale by assets
            diff4A= A/E*N(d1)*sigmaA-sigmaE          # a small number already
            diffNew=abs(diff4E)+abs(diff4A)
            if diffNew<diffOld:
               diffOld=diffNew
               output=(round(A,2),round(sigmaA,4),round(diffNew,5))
    return output
#

def implied_vol_call_min(S,X,T,r,c): 
    from scipy import log,exp,sqrt,stats 
    implied_vol=1.0
    min_value=1000
    for i in range(10000): 
        sigma=0.0001*(i+1)
        d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
        d2 = d1-sigma*sqrt(T)
        c2=S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2) 
        abs_diff=abs(c2-c)
        if abs_diff<min_value: 
            min_value=abs_diff 
            implied_vol=sigma 
            k=i
    return implied_vol


# Step 3: get call option data

def implied_vol_put_min(S,X,T,r,p):
    implied_vol=1.0
    min_value=100.0
    for i in xrange(1,10000): 
        sigma=0.0001*(i+1)
        d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
        d2 = d1-sigma*sqrt(T)
        put=X*exp(-r*T)*stats.norm.cdf(-d2)-S*stats.norm.cdf(-d1) 
        abs_diff=abs(put-p)
        if abs_diff<min_value: 
            min_value=abs_diff 
            implied_vol=sigma 
            k=i
        put_out=put
    print 'k, implied_vol, put, abs_diff' 
    return k,implied_vol, put_out,min_value

def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
#

def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
#

def upCall(s,x,T,r,sigma,nSimulation,barrier):
    import scipy as sp
    import p4f 
    n_steps=100
    dt=T/n_steps 
    inTotal=0 
    outTotal=0
    for j in range(0, nSimulation): 
        sT=s
        inStatus=False 
        outStatus=True
        for i in range(0,int(n_steps)):
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier:
                outStatus=False 
                inStatus=True
        if outStatus==True:
            outTotal+=p4f.bs_call(s,x,T,r,sigma) 
        else:
            inTotal+=p4f.bs_call(s,x,T,r,sigma) 
    return outTotal/nSimulation, inTotal/nSimulation
#

def lookback_min_price_as_strike(s,T,r,sigma,n_simulation): 
    n_steps=100
    dt=T/n_steps
    total=0
    for j in range(n_simulation): 
        min_price=100000.   # a very big number 
        sT=s
        for i in range(int(n_steps)): 
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT<min_price:
                min_price=sT
                #print 'j=',j,'i=',i,'total=',total 
                total+=p4f.bs_call(s,min_price,T,r,sigma)
    return total/n_simulation

#

def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
#

def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100.
    dt=T/n_steps 
    total=0
    for j in range(0, n_simulation): 
        sT=s0
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier:
                out=True 
    if out==False:
       total+=p4f.bs_call(s0,x,T,r,sigma) 
    return total/n_simulation
#

def phaseFirstOrder(self, fromPos, toPos, phi1, degree = True, noChangeOnResonance=False, pivot = 0, scale= "Hz"):
        """This should only be applied to Fourier transformed spectral data
        It will lead to a zero phase shift at the upper end of the spectrum and to a
        phase shift of phi1 at the lower end, linear interpolation inbetween.
        this is the spinsight convention.

        If a pivot is provided the phase will not be changed at the pivot, but the total change accross the entire spectrum will still amount to phi1.
        """
        self.checkToPos(toPos)
        phaseValues = np.linspace(0,phi1,num=len(self.frequency))

        if noChangeOnResonance == True:
            pivot = 0
        elif pivot !=0:
            print("Using pivot for first order phase correction")
            index = self.getIndex(pivot, scale= scale)
            phaseValues = phaseValues - phaseValues[index]
            0
        if degree == True:
            phaseValues = phaseValues*np.pi/180

        self.allFid[toPos] =  [spectrum*np.exp(-1j*phaseValues) for spectrum in self.allFid[fromPos]]

def __MR_W_D_matrix(self,img,labels):
        s = sp.amax(labels)+1
        vect = self.__MR_superpixel_mean_vector(img,labels)

        adj = self.__MR_get_adj_loop(labels)

        W = sp.spatial.distance.squareform(sp.spatial.distance.pdist(vect))

        W = sp.exp(-1*W / self.weight_parameters['delta'])
        W[adj.astype(np.bool)] = 0


        D = sp.zeros((s,s)).astype(float)
        for i in range(s):
            D[i, i] = sp.sum(W[i])

        return W,D

def model(THETA, time, kplanets):
    modelo = 0.0
    if kplanets == 0:
        return 0.0
    for i in range(kplanets):
        As, P, Ac, S, C = THETA[5*i:5*(i+1)]
        A = As ** 2 + Ac ** 2
        ecc = S ** 2 + C ** 2
        w = sp.arccos(C / (ecc ** 0.5))  # longitude of periastron
        phase = sp.arccos(Ac / (A ** 0.5))
        ### test
        if S < 0:
            w = 2 * sp.pi - sp.arccos(C / (ecc ** 0.5))
        if As < 0:
            phase = 2 * sp.pi - sp.arccos(Ac / (A ** 0.5))
        ###
        per = sp.exp(P)
        freq = 2. * sp.pi / per
        M = freq * time + phase
        E = sp.array([MarkleyKESolver().getE(m, ecc) for m in M])
        f = (sp.arctan(((1. + ecc) ** 0.5 / (1. - ecc) ** 0.5) * sp.tan(E / 2.)) * 2.)
        modelo += A * (sp.cos(f + w) + ecc * sp.cos(w))
    return  modelo

def _BlackScholesCall(S, K, T, sigma, r, q):
            d1 = (log(S/K) + (r - q + (sigma**2)/2)*T)/(sigma*sqrt(T))
            d2 = d1 - sigma*sqrt(T)
            return S*exp(-q*T)*norm.cdf(d1) - K*exp(-r*T)*norm.cdf(d2)

def _BlackScholesPut(S, K, T, sigma, r, q):
            d1 = (log(S/K) + (r - q + (sigma**2)/2)*T)/(sigma*sqrt(T))
            d2 = d1 - sigma*sqrt(T)
            return  K*exp(-r*T)*norm.cdf(-d2) - S*exp(-q*T)*norm.cdf(-d1)

def _fprime(self, sigma):
        logSoverK = log(self.S/self.K)
        n12 = ((self.r + sigma**2/2)*self.T)
        numerd1 = logSoverK + n12
        d1 = numerd1/(sigma*sqrt(self.T))
        return self.S*sqrt(self.T)*norm.pdf(d1)*exp(-self.r*self.T)

def kernel(s):
    """
    Constructs a normalized discrete 1D gaussian kernel
    """
    size_grid = int(s*4)
    x= scipy.mgrid[-size_grid:size_grid+1]
    g = scipy.exp(-(x**2/float(s**2)/2.))
    return g / np.sum(g)

def alpha_m(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 0.1*(V+40.0)/(1.0 - sp.exp(-(V+40.0) / 10.0))

def beta_m(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 4.0*sp.exp(-(V+65.0) / 18.0)

def alpha_h(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 0.07*sp.exp(-(V+65.0) / 20.0)

def beta_h(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 1.0/(1.0 + sp.exp(-(V+35.0) / 10.0))

def alpha_n(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 0.01*(V+55.0)/(1.0 - sp.exp(-(V+55.0) / 10.0))

def geom_std(values: t.List[float]) -> float:
    """
    Calculates the geometric standard deviation for the passed values.
    Source: https://en.wikipedia.org/wiki/Geometric_standard_deviation
    """
    import scipy.stats as stats
    import scipy as sp
    gmean = stats.gmean(values)
    return sp.exp(sp.sqrt(sp.sum([sp.log(x / gmean) ** 2 for x in values]) / len(values)))

def asian_vol_adj(atm, time2mat, tau):
    M = (2 * numpy.exp(atm * atm * time2mat) \
         - 2 * numpy.exp(atm * atm * tau) * (1.0 + atm * atm * (time2mat - tau))) / \
        ((atm ** 4) * ((time2mat - tau) ** 2))
    return numpy.sqrt(numpy.log(M) / time2mat)

def BlackSholesFormula(IsCall, S, K, Vol, Texp, Rd, Rf):
    x1 = d1(S, K, Vol, Texp, Rd, Rf )
    x2 = d2(S, K, Vol, Texp, Rd, Rf )
    y = pnorm(x1)
    res = {}

    if IsCall:
        res['Price'] = S * exp(-Rf*Texp)* x1 - K * exp(-Rd*Texp) * x2
        res['Delta'] = x1 * exp(-Rf*Texp)
    else:
        res['Price'] = K * exp(-Rd*Texp) * (1 - x2) - S * exp(-Rf*Texp) * (1 - x1)
        res['Delta'] = (x1  - 1) * exp(-Rf*Texp)
    res['Vega'] = S * sqrt(Texp) * y * exp(-Rf*Texp)
    res['Gamma'] = y * exp(-Rf*Texp)/(S*Vol* sqrt(Texp))
    return res

def KirkApprox(IsCall, F1, F2, Sigma1, Sigma2, Corr, K, Texp, r):
    FA = F1/(F2+K)
    Sigma = sqrt(Sigma1**2 + (Sigma2*F2/(F2+K))**2 - \
                       2*Corr*Sigma1*Sigma2*F2/(F2+K))
    d1 = (numpy.log(FA) + 0.5* Sigma**2 * Texp)/(Sigma*sqrt(Texp))
    d2 = d1 - Sigma*sqrt(Texp)
    x1 = scipy.stats.norm.cdf(d1)
    x2 = scipy.stats.norm.cdf(d2)
    res = {}
    if IsCall:
        res['Price'] = (F2+K)*(FA * x1 - x2) * exp(-r*Texp)
    else:
        res['Price'] = (F2+K)*((1 - x2) - FA*(1 - x1)) * exp(-r*Texp)
    return res

def BSOpt( IsCall, Spot, Strike, Vol, Texp, Rd, Rf ):
    'Standard Black-Scholes European vanilla pricing.'
    if Strike <= 1e-12 * Spot:
        if IsCall:
            return Spot * exp( -Rf * Texp )
        else:
            return 0.

    if IsCall:
        return Spot   * exp( -Rf * Texp ) * cnorm( d1( Spot, Strike, Vol, Texp, Rd, Rf ) ) \
             - Strike * exp( -Rd * Texp ) * cnorm( d2( Spot, Strike, Vol, Texp, Rd, Rf ) )
    else:
        return Strike * exp( -Rd * Texp ) * cnorm( -d2( Spot, Strike, Vol, Texp, Rd, Rf ) ) \
             - Spot   * exp( -Rf * Texp ) * cnorm( -d1( Spot, Strike, Vol, Texp, Rd, Rf ) )

def BSFwd( IsCall, Fwd, Strike, Vol, Texp, ir):
    'Standard Black-Scholes European vanilla pricing.'

    if Strike <= 1e-12 * Fwd:
        if IsCall:
            return Fwd
        else:
            return 0.
    df = exp(-ir * Texp)
    if IsCall:
        return df * (Fwd  * cnorm( fd1( Fwd, Strike, Vol, Texp ) ) \
             - Strike * cnorm( fd2( Fwd, Strike, Vol, Texp ) ))
    else:
        return df * (Strike * cnorm( -fd2( Fwd, Strike, Vol, Texp ) ) \
             - Fwd * cnorm( -fd1( Fwd, Strike, Vol, Texp ) ))

def BSFwdNormal( IsCall, Fwd, Strike, Vol, Texp, ir):
    'Standard Bachelier European vanilla pricing.'
    d = (Fwd-Strike)/Vol/sqrt(Texp)
    p = (Fwd-Strike)  * cnorm( d ) + Vol * sqrt(Texp) * pnorm(d)
    if not IsCall:
        p = p - Fwd + Strike
    return p * exp(-Texp*ir)

def BSDelta( IsCall, Spot, Strike, Vol, Texp, Rd, Rf ):
    'Standard Black-Scholes Delta calculation. Over-currency spot delta.'
    if IsCall:
        return exp( -Rf * Texp ) * cnorm( d1( Spot, Strike, Vol, Texp, Rd, Rf ) )
    else:
        return -exp( -Rf * Texp ) * cnorm( -d1( Spot, Strike, Vol, Texp, Rd, Rf ) )

def BSVega( Spot, Strike, Vol, Texp, Rd, Rf ):
    'Standard Black-Scholes Vega calculation.'

    d = d1( Spot, Strike, Vol, Texp, Rd, Rf )
    return Spot * exp( -Rf * Texp ) * sqrt( Texp / 2. / pi ) * exp( -d * d / 2. )

def BSFwdNormalDelta( IsCall, Fwd, Strike, Vol, Texp, Rd, Rf = 0.0 ):
    d1 = (Fwd - Strike)/Vol/numpy.sqrt(Texp)
    return exp( -Rd * Texp ) * cnorm(d1)

def BSBin( IsCall, Spot, Strike, Vol, Texp, Rd, Rf ):
    'Standard Black-Scholes European binary call/put pricing.'

    Bin = cnorm( d2( Spot, Strike, Vol, Texp, Rd, Rf ) )
    if not IsCall:
        Bin = 1 - Bin
    Bin = Bin * exp( -Rd * Texp )
    return Bin

def StrikeFromDelta( IsCall, Spot, Vol, Texp, Rd, Rf, Delta ):
    '''Calculates the strike of a European vanilla option gives its Black-Scholes Delta.
    It assumes the Delta is an over-ccy spot Delta.'''

    def ArgFunc( Strike ):
        DeltaCalc = BSDelta( IsCall, Spot, Strike, Vol, Texp, Rd, Rf )
        return DeltaCalc - Delta

    LoStrike = Spot * exp( ( Rd - Rf ) * Texp - 4 * Vol * sqrt( Texp ) )
    HiStrike = Spot * exp( ( Rd - Rf ) * Texp + 4 * Vol * sqrt( Texp ) )

    Strike = brenth( ArgFunc, LoStrike, HiStrike )
    return Strike

def OneTouch( IsHigh, IsDelayed, Spot, Strike, Vol, Texp, Rd, Rf ):
    '''Prices a one touch option. IsHigh=True means it knocks up and in; False
    means down and in. IsDelayed=True means it pays at the end; False means it
    pays on hit.'''

    if ( IsHigh and Spot >= Strike ) or ( not IsHigh and Spot <= Strike ):
        if IsDelayed:
            return exp( -Rd * Texp )
        else:
            return 1

    if Vol <= 0 or Texp <= 0: return 0

    Alpha = log( Strike / float( Spot ) )
    Mu    = Rd - Rf - Vol * Vol / 2.

    if IsDelayed:
        if IsHigh:
            Price = exp( -Rd * Texp ) * ( cnorm( ( -Alpha + Mu * Texp ) / Vol / sqrt( Texp ) ) \
                  + exp( 2 * Mu * Alpha / Vol / Vol ) * cnorm( ( -Alpha - Mu * Texp ) / Vol / sqrt( Texp ) ) )
        else:
            Price = exp( -Rd * Texp ) * ( cnorm( (  Alpha - Mu * Texp ) / Vol / sqrt( Texp ) ) \
                  + exp( 2 * Mu * Alpha / Vol / Vol ) * cnorm( (  Alpha + Mu * Texp ) / Vol / sqrt( Texp ) ) )
    else:
        MuHat = sqrt( Mu * Mu + 2 * Rd * Vol * Vol )
        if IsHigh:
            Price = exp( Alpha / Vol / Vol * ( Mu - MuHat ) ) * cnorm( ( -Alpha + MuHat * Texp ) / Vol / sqrt( Texp ) ) \
                  + exp( Alpha / Vol / Vol * ( Mu + MuHat ) ) * cnorm( ( -Alpha - MuHat * Texp ) / Vol / sqrt( Texp ) )
        else:
            Price = exp( Alpha / Vol / Vol * ( Mu + MuHat ) ) * cnorm( (  Alpha + MuHat * Texp ) / Vol / sqrt( Texp ) ) \
                  + exp( Alpha / Vol / Vol * ( Mu - MuHat ) ) * cnorm( (  Alpha - MuHat * Texp ) / Vol / sqrt( Texp ) )

    return Price

def BAWAmOptPricer( IsCall, Fwd, Strike, Vol, Texp, rd, rf ):
    prem = BAWPremium( IsCall, Fwd, Strike, Vol, Texp, rd, rf )
    D    = exp( -rd * Texp)
    Euro = D * BSFwd(IsCall, Fwd, Strike, Vol, Texp)
    return Euro + prem

def LogNormalPaths(mu, cov, fwd, numPaths):
    ''' mu and fwd are 1d lists/arrays (1xn); cov is a 2d scipy.array (nxn); numPaths is int '''
    return (fwd*scipy.exp(numpy.random.multivariate_normal(mu, cov, numPaths) - 0.5*cov.diagonal())).transpose()

def AsianOptTW_Fwd(IsCall, Fwd, strike, RlzAvg, Vol, Texp, AvgPeriod, Rd):
    '''Calculate Asian option price using TurnbullWakeman model.
       This is just for forward options, not for stocks.
       RlzAvg is the realized average price(daily)
       AvgPeriod is the time length of averaging period, usually 22 for the SGX options
    '''

    Fwd = float(Fwd)
    strike = float(strike)
    RlzAvg = float(RlzAvg)

    tau = numpy.max([0, Texp - AvgPeriod])
    if AvgPeriod == 0:
        volA = Vol
    else:
        volA = asian_vol_adj(Vol, Texp, tau)

    X = numpy.copy(strike)
    if AvgPeriod > Texp:
        X = X * (AvgPeriod / Texp) - RlzAvg * (AvgPeriod - Texp) / Texp

    if X < 0:
        if IsCall:
            return (RlzAvg * (AvgPeriod - Texp) / AvgPeriod + Fwd * Texp / AvgPeriod - X) * exp(-Rd * Texp)
        else:
            return 0
    else:
        price = BSFwd(IsCall, Fwd, X, volA, Texp, Rd)
        if AvgPeriod > Texp:
            return price * Texp / AvgPeriod
        else:
            return price

def AsianFwdDelta(IsCall, Fwd, strike, RlzAvg, Vol, Texp, AvgPeriod, Rd):
    Fwd = float(Fwd)
    strike = float(strike)
    RlzAvg = float(RlzAvg)

    if AvgPeriod > Texp:
        x = strike * (AvgPeriod / Texp) - RlzAvg * (AvgPeriod - Texp) / Texp
    else:
        x = strike
    tau = numpy.max([0, Texp - AvgPeriod])

    if AvgPeriod > 0:
        volA = asian_vol_adj(Vol, Texp, tau)
    else:
        volA = Vol

    if x < 0:
        if IsCall:
            return exp(-Rd * Texp) * Texp / AvgPeriod
        else:
            return 0
    else:
        if Texp < AvgPeriod:
            multi = Texp / AvgPeriod
        else:
            multi = 1.0
        Asiand1 = (log(Fwd / x) + (volA * volA * 0.5) * Texp) / (volA * sqrt(Texp))
        if IsCall:
            return multi * exp(-Rd * Texp) * cnorm(Asiand1)
        else:
            return multi * exp(-Rd * Texp) * (cnorm(Asiand1) - 1.0)

def AsianFwdGamma(Fwd, strike, RlzAvg, Vol, Texp, AvgPeriod, Rd):
    Fwd = float(Fwd)
    strike = float(strike)
    RlzAvg = float(RlzAvg)

    if AvgPeriod > Texp:
        x = strike * (AvgPeriod / Texp) - RlzAvg * (AvgPeriod - Texp) / Texp
    else:
        x = strike
    tau = numpy.max([0, Texp - AvgPeriod])

    if AvgPeriod > 0:
        volA = asian_vol_adj(Vol, Texp, tau)
    else:
        volA = Vol

    if x < 0:
        return 0
    else:
        if Texp < AvgPeriod:
            multi = Texp / AvgPeriod
        else:
            multi = 1.0

        Asiand1 = (log(Fwd / x) + (volA * volA * 0.5) * Texp) / (volA * sqrt(Texp))
        ND = exp(-(Asiand1 * Asiand1 * 0.5)) / sqrt(2 * pi)
        return multi * exp(- Rd * Texp) * ND / (Fwd * volA * sqrt(Texp))

def AsianFwdVega(Fwd, strike, RlzAvg, Vol, Texp, AvgPeriod, Rd):
    Fwd = float(Fwd)
    strike = float(strike)
    RlzAvg = float(RlzAvg)

    if AvgPeriod > Texp:
        x = strike * (AvgPeriod / Texp) - RlzAvg * (AvgPeriod - Texp) / Texp
    else:
        x = strike
    tau = numpy.max([0, Texp - AvgPeriod])

    if AvgPeriod > 0:
        M = (2.0 * exp(Vol * Vol * Texp) - 2.0 * exp(Vol * Vol * tau) * (1.0 + Vol * Vol * (Texp - tau))) / \
            ((Vol ** 4) * ((Texp - tau) ** 2))
        volA = sqrt(log(M) / Texp)
    else:
        volA = Vol

    if x < 0:
        return 0
    else:
        if Texp < AvgPeriod:
            multi = Texp / AvgPeriod
        else:
            multi = 1.0
        if AvgPeriod > 0:
            dM = 4.0 * (exp(Vol * Vol * Texp) * Texp * Vol - exp(Vol * Vol * tau) \
                        * ((Vol ** 3) * tau * (Texp - tau) + Vol * Texp)) / \
                 ((Vol ** 4) * (Texp - tau) * (Texp - tau)) - \
                 8.0 * (exp(Vol * Vol * Texp) - exp(Vol * Vol * tau) * (1.0 + Vol * Vol * (Texp - tau))) / \
                 ((Vol ** 5) * (Texp - tau) * (Texp - tau))

            dvA = 1.0 / (2.0 * volA) / Texp / M * dM
        else:
            dvA = 1.0
        Asiand1 = (log(Fwd / x) + volA * volA * 0.5 * Texp) / (volA * sqrt(Texp))
        ND = exp(-(Asiand1 * Asiand1 * 0.5)) / sqrt(2 * pi)
        return multi * Fwd * exp(-Rd * Texp) * ND * sqrt(Texp) * dvA * 0.01

def devi(yy, eta):
    deveta = yy*eta - scipy.exp(eta)
    devy = yy*scipy.log(yy) - yy
    devy[yy == 0] = 0
    result = 2*(devy - deveta)
    return(result)

def devi(yy, eta):
    deveta = yy*eta - scipy.exp(eta)
    devy = yy*scipy.log(yy) - yy
    devy[yy == 0] = 0
    result = 2*(devy - deveta)
    return(result)

def devi(yy, eta):
    deveta = yy*eta - scipy.exp(eta)
    devy = yy*scipy.log(yy) - yy
    devy[yy == 0] = 0
    result = 2*(devy - deveta)
    return(result)

def kde(data, N=None, MIN=None, MAX=None):

    # Parameters to set up the mesh on which to calculate
    N = 2**14 if N is None else int(2**sci.ceil(sci.log2(N)))
    if MIN is None or MAX is None:
        minimum = min(data)
        maximum = max(data)
        Range = maximum - minimum
        MIN = minimum - Range/10 if MIN is None else MIN
        MAX = maximum + Range/10 if MAX is None else MAX

    # Range of the data
    R = MAX-MIN

    # Histogram the data to get a crude first approximation of the density
    M = len(data)
    DataHist, bins = sci.histogram(data, bins=N, range=(MIN,MAX))
    DataHist = DataHist/M
    DCTData = scipy.fftpack.dct(DataHist, norm=None)

    I = [iN*iN for iN in xrange(1, N)]
    SqDCTData = (DCTData[1:]/2)**2

    # The fixed point calculation finds the bandwidth = t_star
    guess = 0.1
    try:
        t_star = scipy.optimize.brentq(fixed_point, 0, guess, 
                                       args=(M, I, SqDCTData))
    except ValueError:
        print 'Oops!'
        return None

    # Smooth the DCTransformed data using t_star
    SmDCTData = DCTData*sci.exp(-sci.arange(N)**2*sci.pi**2*t_star/2)
    # Inverse DCT to get density
    density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R
    mesh = [(bins[i]+bins[i+1])/2 for i in xrange(N)]
    bandwidth = sci.sqrt(t_star)*R

    density = density/sci.trapz(density, mesh)
    return bandwidth, mesh, density

def fixed_point(t, M, I, a2):
    l=7
    I = sci.float128(I)
    M = sci.float128(M)
    a2 = sci.float128(a2)
    f = 2*sci.pi**(2*l)*sci.sum(I**l*a2*sci.exp(-I*sci.pi**2*t))
    for s in range(l, 1, -1):
        K0 = sci.prod(xrange(1, 2*s, 2))/sci.sqrt(2*sci.pi)
        const = (1 + (1/2)**(s + 1/2))/3
        time=(2*const*K0/M/f)**(2/(3+2*s))
        f=2*sci.pi**(2*s)*sci.sum(I**s*a2*sci.exp(-I*sci.pi**2*time))
    return t-(2*M*sci.sqrt(sci.pi)*f)**(-2/5)

def computeOpenMaxProbability(openmax_fc8, openmax_score_u):
    """ Convert the scores in probability value using openmax

    Input:
    ---------------
    openmax_fc8 : modified FC8 layer from Weibull based computation
    openmax_score_u : degree

    Output:
    ---------------
    modified_scores : probability values modified using OpenMax framework,
    by incorporating degree of uncertainity/openness for a given class

    """
    prob_scores, prob_unknowns = [], []
    for channel in range(NCHANNELS):
        channel_scores, channel_unknowns = [], []
        for category in range(NCLASSES):
            channel_scores += [sp.exp(openmax_fc8[channel, category])]

        total_denominator = sp.sum(sp.exp(openmax_fc8[channel, :])) + sp.exp(sp.sum(openmax_score_u[channel, :]))
        prob_scores += [channel_scores/total_denominator ]
        prob_unknowns += [sp.exp(sp.sum(openmax_score_u[channel, :]))/total_denominator]

    prob_scores = sp.asarray(prob_scores)
    prob_unknowns = sp.asarray(prob_unknowns)

    scores = sp.mean(prob_scores, axis = 0)
    unknowns = sp.mean(prob_unknowns, axis=0)
    modified_scores =  scores.tolist() + [unknowns]
    assert len(modified_scores) == 1001
    return modified_scores

#---------------------------------------------------------------------------------

def __init__(self, data, mleValue, fitParameters=True, mean=None, sigma=None):
        super(normalModel, self).__init__(data)




        self.MLE = mleValue

        if(None in [mean, sigma]):
            fitParameters=True
        if(fitParameters):
            mean = Mean(self.getDataSet())      
            sigma = standardDeviation(self.getDataSet())
            try:

                def normDist(x, x0, sigma):
                    output = 1.0/sqrt(2*np.pi*(sigma**2))
                    output *= exp(-0.5*((x - x0)/sigma)**2)
                    return output

                self.n, self.bins, patches = plt.hist(self.getDataSet(), self.getDatasetSize()/10, normed=1, facecolor='blue', alpha = 0.55)
                popt,pcov = curve_fit(normDist,self.bins[:-1], self.n, p0=[mean, sigma])
                ##plt.plot(bins[:-1], gaus(bins[:-1],*popt),'c-',label="Gaussian Curve with params\na=%f\nx0=%f\nsigma=%f" % (popt[0], popt[1], popt[2]), alpha=0.5)
                print "Fitted gaussian curve to data with params x0 %f, sigma %f" % (popt[0], popt[1])
                self.x0 = popt[0]
                self.sigma = popt[1]

                self.fitted = True
            except RuntimeError:
                print "Unable to fit data to normal curve, runtime error"
                raise
            except Warning:
                raise RuntimeError
        else:
            self.x0 = mean
            self.sigma = sigma

def getModelpdf(self, x):
        output = 1.0/math.sqrt(2*np.pi*(self.getSigmaValue()**2))
        output *= math.exp(-0.5*((x - self.getx0Value())/self.getSigmaValue())**2)
        return output