brownian bridge and quasi-MC

[sdementen]

first things first, great package, tx! I have just tried some examples but look very promising.

what is the control one can have on the generation of random numbers while generating tracjectories ?
I would like to be able to use a quasi-MC (sobol) RNG combined with some brownian bridge transformation. Is there a hook I can use to provide a function that generates the random numbers ?

[RadostW]

Thanks! 😊

First of all, you can already get Brownian bridge samples with the current backend. It's not immediately obvious but Brownian bridge is actually an Ito process itself

(cf https://math.stackexchange.com/questions/1369834/brownian-bridge-sde )

If you want to drive your process with Brownian bridge then with some algebra you can express your SDE against dB as an SDE against dW (with B bridge and W wiener processes respectively).

Take a look at the code below:

import pychastic
import jax.numpy as jnp

a = 1
b = -1
tmax = 1

problem = pychastic.sde_problem.SDEProblem(
    lambda x: jnp.array([1, (b - x[1]) / (1 - x[0])]),
    lambda x: jnp.array([[0], [1]]),
    x0=jnp.array([0, a]),
    tmax=tmax,
)
solver = pychastic.sde_solver.SDESolver(dt=0.001)
trajectory = solver.solve_many(problem,n_trajectories = 10)

import matplotlib.pyplot as plt

for t in trajectory["solution_values"]:
    plt.plot(trajectory["time_values"][0,:], t[:,1])
plt.show()

If you really need to control the RNG then it's controlled by
get_wiener_integrals function from:
"https://github.com/RadostW/stochastic/blob/main/pychastic/vectorized_I_generation.py"
To support higher order schemes this function is rather complex because we need approximations of multiple integrals.

I'm not familiar with "quasi-MC (sobol)" what does that mean? Could you provide more detail on the problem you're trying to solve?

[sdementen]

Great for BB, I'll have a look. For sobol and other quasi random numbers, the goal is to get 1/N convergence on monte Carlo instead of 1/sqrt(N). I have attached a doc where the technique is nicely described.

…

On Thu, 5 Oct 2023, 16:39 RadostW, ***@***.***> wrote: Thanks! 😊 First of all, you can already get Brownian bridge samples with the current backend. It's not immediately obvious but Brownian bridge is actually an Ito process itself (cf https://math.stackexchange.com/questions/1369834/brownian-bridge-sde ) If you want to drive your process with Brownian bridge then with some algebra you can express your SDE against dB as an SDE against dW (with B bridge and W wiener processes respectively). Take a look at the code below: import pychasticimport jax.numpy as jnp a = 1b = -1tmax = 1 problem = pychastic.sde_problem.SDEProblem( lambda x: jnp.array([1, (b - x[1]) / (1 - x[0])]), lambda x: jnp.array([[0], [1]]), x0=jnp.array([0, a]), tmax=tmax, )solver = pychastic.sde_solver.SDESolver(dt=0.001)trajectory = solver.solve_many(problem,n_trajectories = 10) import matplotlib.pyplot as plt for t in trajectory["solution_values"]: plt.plot(trajectory["time_values"][0,:], t[:,1])plt.show() If you *really* need to control the RNG then it's controlled by get_wiener_integrals function from: " https://github.com/RadostW/stochastic/blob/main/pychastic/vectorized_I_generation.py " To support higher order schemes this function is rather complex because we need approximations of multiple integrals. I'm not familiar with "quasi-MC (sobol)" what does that mean? Could you provide more detail on the problem you're trying to solve? — Reply to this email directly, view it on GitHub <#10 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAJ6S5TUZTRESJ5OSHBAGQDX53BBDAVCNFSM6AAAAAA5UHKZ7CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONBZGA2DIMRYHE> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[RadostW]

I think attachments don't work properly inside GitHub threads.

The Sobol sequence idea sounds very interesting - it seems scipy has the initialization numbers done correctly up to dimension 21201. I haven't looked very deep but it seems that we could just draw uniformly on a hypercube and then transform with normal CDF.

Perhaps an integration scheme "sobol-euler" could be implemented relatively easily and exposed to the users similarly to other solvers?

Inside the package there are already some correctness tests so debugging should go smoothly.

[sdementen]

The link to the PDF https://personalpages.manchester.ac.uk/staff/paul.johnson-2/resources/math60082/notes-math60082-3.pdf The sobol sequence is best combined with Brownian bridge to accelerate convergence of Montecarlo method, so ideally, the two must be combined

…

On Thu, 5 Oct 2023, 21:47 RadostW, ***@***.***> wrote: I think attachments don't work properly inside GitHub threads. The Sobol sequence idea sounds very interesting - it seems scipy has the initialization numbers done correctly up to dimension 21201. I haven't looked very deep but it seems that we could just draw uniformly on a hypercube and then transform with normal CDF. Perhaps an integration scheme "sobol-euler" could be implemented relatively easily and exposed to the users similarly to other solvers? Inside the package there are already some correctness tests so debugging should go smoothly. — Reply to this email directly, view it on GitHub <#10 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAJ6S5X7JKTUIG3YGKROOA3X54FD7AVCNFSM6AAAAAA5UHKZ7CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONBZGU2DGMZWGY> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[sdementen]

Some other papers/slides on the topic I found deeply interesting and well written: https://people.cs.kuleuven.be/~dirk.nuyens/taiwan/Talk_Kritzer.pdf https://www.ynufe.edu.cn/pub/tsxy/jhzhao/teach/CS/notes/lec5.pdf https://people.maths.ox.ac.uk/gilesm/mc/module_6/QMC.pdf

…

On Thu, 5 Oct 2023, 22:46 Sébastien de Menten, ***@***.***> wrote: The link to the PDF https://personalpages.manchester.ac.uk/staff/paul.johnson-2/resources/math60082/notes-math60082-3.pdf The sobol sequence is best combined with Brownian bridge to accelerate convergence of Montecarlo method, so ideally, the two must be combined On Thu, 5 Oct 2023, 21:47 RadostW, ***@***.***> wrote: > I think attachments don't work properly inside GitHub threads. > > The Sobol sequence idea sounds very interesting - it seems scipy has the > initialization numbers done correctly up to dimension 21201. I haven't > looked very deep but it seems that we could just draw uniformly on a > hypercube and then transform with normal CDF. > > Perhaps an integration scheme "sobol-euler" could be implemented > relatively easily and exposed to the users similarly to other solvers? > > Inside the package there are already some correctness tests so debugging > should go smoothly. > > — > Reply to this email directly, view it on GitHub > <#10 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AAJ6S5X7JKTUIG3YGKROOA3X54FD7AVCNFSM6AAAAAA5UHKZ7CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONBZGU2DGMZWGY> > . > You are receiving this because you authored the thread.Message ID: > ***@***.***> >

[sdementen]

I have python code that builds up the matrix L in QMC.pdf p 25, 26 and 29 to use classic MC, PCA or BB