#!/usr/bin/env python3
"``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
)
