How to simulate a stochastic process until all elements of a path are positive?

Question

I want to simulate a path for certain a stochastic process in continuous time using Python. I wrote the following function to simulate it (np refers to numpy of course):

def simulate_V(v0, psi, theta, dt, xsi, sample, diff = False):

    _, dB = independent_brownian(dt, 1, sample, diff = True)
    path_lenght = len(dB) + 1
    V = np.zeros(path_lenght)
    dV = np.zeros(len(dB))
    V[0] = v0
    for i in range(len(dV)):
        dV[i] = psi * (theta - V[i]) * dt + xsi * np.sqrt(V[i]) * dB[i]
        V[i+1] = V[i] + dV[i]
    
    if diff == True:
        return(V, dV)
    else:
        return(V)

The independent_brownian function just creates a path for a standard Brownian motion. For completeness, here it is:

def independent_brownian(dt, n, sample, diff=False):
    '''
    Creates paths of independent Brownian motions. Returns increments if diff == True.
    --------
    dt -> variance of increments\\
    n -> how many paths\\
    sample -> length of sampled path\\
    diff -> return increments or not?\\

    Returns a (n, sample)-shaped array with the paths
    '''
    increments = np.sqrt(dt) * np.random.randn(n, sample)
    paths = np.cumsum(increments, axis=1)
    b = np.zeros((n, sample + 1))
    b[:, 1:] = paths

    if n==1:
        b = b.flatten()
        increments = increments.flatten()

    if diff == True:
        return(b, increments)
    else:
        return(b)

It happens that the math behind my model implies that the process $V_t$, which is represented by its discretization V in the code above, must be positive. However, it might be the case where the array dB contains negative elements which are big in absolute value. I would like to automate the following "choice" procedure:

1. Try to follow the loop described inside simulate_V;
2. At some point, if V[i] falls below zero, interrupt the process, sample another sequence dB and start over;
3. Stop when all elements of V are positive;

What is a good way to automate this procedure? Right now, I understand that, as it is, if V[i] goes below zero, I will get a nan in numpy but it will throw no error or stop the process.

Answer 1

For the sake of generality, and since I am not very familiar with Brownian Motion, I will implement a general mechanism for the situation in which we have a variable dB, which we use to create another variable V, and we need to repeat this generation of V until a certain condition is held.

import numpy as np
dB = np.random.normal(size=3) # initializing dB
found_sufficient_v_flag = False # we define a flag that would tell us when to stop
while not found_sufficient_v_flag: # while this flag is False, we continue searching
  v = np.zeros(3) # initializing v
  for i in range(3): # for loop
    v[i] = dB[i] + np.random.normal() # we change the v[i] element
    if v[i] < 0: # if it is negative, there is no point in continuing, so we break the loop
      dB = np.random.normal(size=3) # but we need to instantiate a different dB
      break 
  if np.all(v > 0): # we arrive here if either v[i] < 0, or all v[i] > 0. in the latter case we stop
    found_sufficient_v_flag = True
print(dB)
print(v)

This provided me with the following output:

[2.27582634 0.77008881 0.28388536]
[2.55101104 3.10944337 0.55829105]

And you can see that the condition on V does hold.

If any other clarification is needed, let me know.