Source code for naginterfaces.​library.​examples.​rand.​bb_inc_ex

#!/usr/bin/env python
"``naginterface.library.rand.bb_inc`` Python Example."

# NAG Copyright 2017-2019.

# pylint: disable=invalid-name,too-many-locals

from math import sqrt

import numpy as np

from naginterfaces.library import rand

[docs]def main(): """ Example for :func:`naginterfaces.library.rand.bb_inc`. Backs out increments from sample paths generated by a Brownian bridge algorithm. >>> main() naginterfaces.library.rand.bb_inc Python Example Results. Compute numerical solutions to an SDE driven by a Wiener process. Euler--Maruyama solution for Geometric Brownian motion: Path 1 Path 2 Path 3 1 0.9668 1.0367 0.9992 2 0.9717 1.0254 1.0077 3 0.9954 1.0333 1.0098 4 0.9486 1.0226 0.9911 5 0.9270 1.0113 1.0630 6 0.8997 1.0127 1.0164 7 0.8955 1.0138 1.0771 8 0.8953 0.9953 1.0691 9 0.8489 1.0462 1.0484 10 0.8449 1.0592 1.0429 11 0.8158 1.0233 1.0625 12 0.7997 1.0384 1.0729 13 0.8025 1.0138 1.0725 14 0.8187 1.0499 1.0554 15 0.8270 1.0459 1.0529 16 0.7914 1.0294 1.0783 17 0.8076 1.0224 1.0943 18 0.8208 1.0359 1.0773 19 0.8190 1.0326 1.0857 20 0.8217 1.0326 1.1095 21 0.8084 0.9695 1.1389 Analytic solution at final time step: Path 1 Path 2 Path 3 0.8079 0.9685 1.1389 """ print('naginterfaces.library.rand.bb_inc Python Example Results.') print('Compute numerical solutions to an SDE driven by a Wiener process.') # The time discretization on which to solve the SDE: t0 = 0. tend = 1. ntimes = 20 # The interior time points: intime = np.linspace(1./(ntimes+1), tend, num=ntimes, endpoint=False) # Create the Brownian bridge construction order; # order 3, do not move and times to the front: times = rand.bb_make_bridge_order(t0, tend, intime, bgord=3) # Initialize the Brownian bridge increments generator: bb_comm = rand.bb_inc_init(t0, tend, times) # The configuration for computing the scaled increments of the Wiener # process: npaths = 3 d = 1 # The dimension of each sample path. a = 0 # Using a free Wiener process. # Terminal/starting differences not required for a free Wiener process: diff = [0.] # The lower-triangular Cholesky factorization for the covariance # matrix of the Wiener process: c = [[1.]] # Generate the Normal data for constructing the sample paths. # Initialize the base PRNG to a repeatable sequence; # use L'Ecuyer MRG32k3a: genid = 6 statecomm = rand.init_repeat(genid, seed=[1023401]) # Initialize the scrambled QRNG (scrambled quasi-random sequences # counteract the bias some applications experience when using # quasi-random sequences); # use Sobol: genid = 1 # Use Faure--Tezuka scrambling: stype = 2 # The number of quasi-random points required: idim = d*(ntimes + 1 - a) # Skip no terms on initialization: iskip = 0 # Scramble all digits: nsdigi = 0 sqrng_comm = rand.quasi_init_scrambled( genid, stype, idim, iskip, nsdigi, statecomm, ) # Generate npaths quasi-random points from N(0,1): xmean = [0.]*idim std = [1.]*idim z = rand.quasi_normal(xmean, std, npaths, sqrng_comm) # Get the scaled increments of the Wiener process: sinc = rand.bb_inc(npaths, diff, z, c, bb_comm, a=a) # Do the Euler--Maruyama time stepping for the SDE: s0 = 1. r = 0.05 sigma = 0.12 st = np.empty((npaths, ntimes+1)) st[:, 0] = s0 + (intime[0]-t0)*(r*s0+sigma*s0*sinc.b[:, 0]) for i in range(1, ntimes): st[:, i] = ( st[:, i-1] + (intime[i]-intime[i-1])* (r*st[:, i-1]+sigma*st[:, i-1]*sinc.b[:, i]) ) st[:, ntimes] = ( st[:, ntimes-1] + (tend-intime[ntimes-1])* (r*st[:, ntimes-1]+sigma*st[:, ntimes-1]*sinc.b[:, ntimes]) ) # The analytic solution is # ST = S0*exp( (r-sigma**2/2)T + sigma WT ) analytic = s0*np.exp((r-sigma**2/2.)*tend+sigma*sqrt(tend-t0)*sinc.z[:, 0]) print('Euler--Maruyama solution for Geometric Brownian motion:') path_head_str = ( ' '*4 + ''.join([ ' '*5 + 'Path {:d}'.format(p+1) for p in range(npaths) ]) ) print(path_head_str) for i in range(ntimes + 1): print( ' {:2d}'.format(i+1) + ' ' + ''.join([' {:10.4f}']*npaths).format(*st[:, i]) ) print('Analytic solution at final time step:') print(path_head_str) print(' '*4 + ''.join([' {:10.4f}']*npaths).format(*analytic))
if __name__ == '__main__': import doctest import sys sys.exit( doctest.testmod( None, verbose=True, report=False, optionflags=doctest.REPORT_NDIFF, ).failed )