-
Notifications
You must be signed in to change notification settings - Fork 1
/
brownian.py
105 lines (91 loc) · 2.9 KB
/
brownian.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
from __future__ import division
import math
import scipy.stats
def Phi(x):
"""
Cumulative standard normal distribution.
"""
return scipy.stats.norm.cdf(x)
def U(n,gamma,A,y):
"""
This is a primitive function of e^(gamma y)sin ((n pi y)/A),
multiplied by 2/A*exp(-gamma*y).
"""
return (2*A*gamma*math.sin(math.pi*n*y/A) - 2*math.pi*n*math.cos (math.pi*n*y/A))/(A**2*gamma**2 + math.pi**2 * n**2)
class Brownian:
def __init__(self,a=-1.0,b=1.0,mu=0.0,sigma=0.005):
self.a=a
self.b=b
self.mu=mu
self.sigma=sigma
self.sigma2=sigma**2
gamma=self.mu/self.sigma2
def outcome_cdf(self,T=None,y=None):
# in case of slow convergence use Siegmund approximation.
sigma2=self.sigma2
mu=self.mu
gamma=mu/sigma2
A=self.b-self.a
if sigma2*T/A**2<1e-2 or abs(gamma*A)>15:
ret=self.outcome_cdf_alt2(T,y)
else:
ret=self.outcome_cdf_alt1(T,y)
assert -1e-3 <= ret <= 1+1e-3
return ret
def outcome_cdf_alt1(self,T=None,y=None):
"""
Computes the probability that the particle passes to the
right of (T,y), the time axis being vertically oriented.
This may give a numerical exception if math.pi**2*sigma2*T/(2*A**2)
is small.
"""
mu=self.mu
sigma2=self.sigma2
A=self.b-self.a
x=0-self.a
y=y-self.a
gamma=mu/sigma2
n=1
s=0.0
lambda_1=((math.pi/A)**2)*sigma2/2+(mu**2/sigma2)/2
t0=math.exp(-lambda_1*T-x*gamma+y*gamma)
while True:
lambda_n=((n*math.pi/A)**2)*sigma2/2+(mu**2/sigma2)/2
t1=math.exp(-(lambda_n-lambda_1)*T)
t3=U(n,gamma,A,y)
t4=math.sin(n*math.pi*x/A)
s+=t1*t3*t4
if abs(t0*t1*t3)<=1e-9:
break
n+=1
if gamma*A>30: # avoid numerical overflow
pre=math.exp(-2*gamma*x)
elif abs(gamma*A)<1e-8: # avoid division by zero
pre=(A-x)/A
else:
pre=(1-math.exp(2*gamma*(A-x)))/(1-math.exp(2*gamma*A))
return pre+t0*s
def outcome_cdf_alt2(self,T=None,y=None):
"""
Siegmund's approximation. We use it as backup if our
exact formula converges too slowly. To make the evaluation
robust we use the asymptotic development of Phi.
"""
denom=math.sqrt(T*self.sigma2)
offset=self.mu*T
gamma=self.mu/self.sigma2
a=self.a
b=self.b
z=(y-offset)/denom
za=(-y+offset+2*a)/denom
zb=(y-offset-2*b)/denom
t1=Phi(z)
if gamma*a>=5:
t2=-math.exp(-za**2/2+2*gamma*a)/math.sqrt(2*math.pi)*(1/za-1/za**3)
else:
t2=math.exp(2*gamma*a)*Phi(za)
if gamma*b>=5:
t3=-math.exp(-zb**2/2+2*gamma*b)/math.sqrt(2*math.pi)*(1/zb-1/zb**3)
else:
t3=math.exp(2*gamma*b)*Phi(zb)
return t1+t2-t3