#!/usr/bin/python -Wall # ================================================================ # John Kerl # kerl at math dot arizona dot edu # 2008-05-12 # ================================================================ from __future__ import division # 1/2 = 0.5, not 0. import random from math import * import copy import sys import re # ---------------------------------------------------------------- def usage(): print >> sys.stderr, "Usage: %s [options]" % (sys.argv[0]) print >> sys.stderr, "Options:" print >> sys.stderr, " T=[number]" print >> sys.stderr, " nt=[number]" print >> sys.stderr, " X0=[number]" print >> sys.stderr, " m=[number]" print >> sys.stderr, " sigma=[number]" print >> sys.stderr, " nX=[number]" sys.exit(1) # ---------------------------------------------------------------- # Mean-reverting Ornstein-Uhlenbeck process: # SDE: dX = (m-X) dt + sigma db # Mean of solution: X0 exp(-t) + m(1 - exp(-t)). T = 10 nt = 1000 X0 = 1.7 m = 2 sigma = 0.05 nX = 6 argc = len(sys.argv) for argi in range(1, argc): if re.match(r'^T=', sys.argv[argi]): T = float(sys.argv[argi][2:]) elif re.match(r'^nt=', sys.argv[argi]): nt = int (sys.argv[argi][3:]) elif re.match(r'^X0=', sys.argv[argi]): X0 = float(sys.argv[argi][3:]) elif re.match(r'^m=', sys.argv[argi]): m = float(sys.argv[argi][2:]) elif re.match(r'^sigma=', sys.argv[argi]): sigma = float(sys.argv[argi][6:]) elif re.match(r'^nX=', sys.argv[argi]): nX = int (sys.argv[argi][3:]) else: usage() t = 0 dt = T/nt X = [X0] * nX Bt = [0.0] * nX sqrtdt = sqrt(dt) while (t <= T): mean = X0 * exp(-t) + m*(1 - exp(-t)) stddev = sqrt(sigma**2/2 * (1 - exp(-2*t))) print "%11.7f %11.7f %11.7f %11.7f %11.7f %11.7f" % \ (t, mean+2*stddev, mean+stddev, mean, mean-stddev, mean-2*stddev), for k in range(0, nX): print "%11.7f" % (X[k]), print for k in range(0, nX): dB = random.gauss(0, sqrtdt) X[k] += (m - X[k]) * dt + sigma * dB Bt[k] += dB t += dt