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 ) )
python类exp()的实例源码
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
