#!/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, " c=[number]" print >> sys.stderr, " alpha=[number,number,...]" print >> sys.stderr, " nX=[number]" sys.exit(1) # ---------------------------------------------------------------- # Process: X_t = exp(c t + sum_{j=1}^n alpha_j btj) # # = exp(c t) * prod_{j=1}^n exp(alpha_j btj) # # There are nX realizations plotted here. Each of them has n Brownian # motions attached: Bt[k][j] for k from 0 to nX-1 and j from 0 to n-1. T = 2 nt = 2000 c = 0.4 alpha = [0.1, -0.3, 0.2] nX = 20 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'^c=', sys.argv[argi]): c = float(sys.argv[argi][2:]) elif re.match(r'^alpha=', sys.argv[argi]): alpha = [ float(s) for s in sys.argv[argi][6:].split(',') ] elif re.match(r'^nX=', sys.argv[argi]): nX = int (sys.argv[argi][3:]) else: usage() t = 0 dt = T/nt X0 = 1 n = len(alpha) X = [X0] * nX Bt = [] for k in range(0, nX): Bt.append(copy.copy([0.0] * n)) sqrtdt = sqrt(dt) while (t <= T): for k in range(0, nX): X[k] = exp(c*t) for j in range(0, n): X[k] *= exp(alpha[j] * Bt[k][j]) print "%11.7f" % (t), for k in range(0, nX): print "%11.7f" % (X[k]), print for k in range(0, nX): for j in range(0, n): dB = random.gauss(0, sqrtdt) Bt[k][j] += dB t += dt