Guest User

Untitled

a guest
Jun 23rd, 2018
86
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.82 KB | None | 0 0
  1. # Use Euler-Maruyama to generate realizations of Brownian Bridge
  2. # stochastic differential equations
  3. #
  4. # Euler-Maruyama method for Stochastic Differential Equation
  5. # dy(t)=f(t,y)dt+g(t,y)dB
  6. # Input: functions f and g, initial y0, arrays t and B
  7. # t contains times t[i], i=0 to n at which solution desired
  8. # B contains B[i], i=0 to n basic Brownian motion at these times
  9. # Output: array w with approximate solution w[i], i=0 to n
  10. from numpy import *
  11. import numpy as np
  12. from random import *
  13.  
  14. def eulermaruyama(f,g,y0,t,B):
  15. n = shape(t)[0]-1; w = empty([n+1]); w[0] = y0
  16. for i in range(n):
  17. Deltati=t[i+1]-t[i]; DeltaBi = B[i+1]-B[i]
  18. w[i+1]=w[i]+f(t[i],w[i])*Deltati+g(t[i],w[i])*DeltaBi
  19. return w
  20.  
  21. #
  22.  
  23.  
  24. (I THINK this following might be what was use, too.)
  25.  
  26. #
  27. a=1.; b=3.; y0=1. # problem parameters
  28. N=1000 # N large is number of subintervals for basic Brownian motion
  29. t=linspace(a,b,N+1,endpoint=True); # mesh both for basic Brownian motion and for E-M method approx.
  30. nrealizations=4
  31. w=empty([nrealizations,N+1]); B=empty([nrealizations,N+1])
  32. colors=array(['r','g','b','k'])
  33. for i in range(nrealizations):
  34. z=np.random.randn(N) # choice of random variables zz
  35. B[i,:]=basic_brownian(t,z) # determines basic Brownian
  36. w[i,:]=eulermaruyama(f,g,y0,t,B[i,:]) # Euler-Maruyama approximate solution at these times
  37.  
  38. import matplotlib.pyplot as plt
  39. for i in range(nrealizations):
  40. plt.plot(t,B[i,:],colors[i]+'-',lw=0.5) # lw specifies line width, 1.0 is default
  41. plt.plot(t,w[i,:],colors[i]+'-',lw=1.5)
  42. plt.xlabel('t')
  43. plt.title('Figure 9.16'+("%2d" % nrealizations)+' solutions (bold) to the Brownian bridge equation')
  44. plt.savefig('eulermaruyama_brownianbridge.png')
  45. plt.show()
Add Comment
Please, Sign In to add comment