Simulating an elementary stochastic process in Python

Question

I'm trying to simulate a simple stochastic process in Python, but with no success. The process is the following:

x(t + δt) = r(t) * x(t)

where r(t) is a Bernoulli random variable that can assume the values 1.5 or 0.6.

I've tried the following:

n = 10
r = np.zeros( (1,n))

for i in range(0, n, 1):
    if r[1,i] == r[1,0]:
        r[1,i] = 1
    else:
        B = bernoulli.rvs(0.5, size=1)
        if B == 0:
            r[1,i] = r[1,i-1] * 0.6
        else:
            r[1,i] = r[1,i-1] * 1.5

Can you explain why this is wrong and a possible solution?

The code regarding the Bernoulli random variable is confusing, please encapsulate it into a function `get_bernoulli` that returns either 0.6 or 1.5 — peer, Feb 14 '21 at 18:01
Python is zero indexed, so perhaps the first index to `r` should be 0? — Alex, Feb 14 '21 at 20:20

score 0 · Answer 1 · answered Sep 24 '21 at 16:03

So , first thing is that the SDE should be perceived over time, so you also need to consider the discretization rather than just giving the number of steps through n .

Essentially, what you are asking is just a simple random walk with a Bernoulli random variable taking on the values 0.5 and 1.6 instead of a Gaussian (standard normal) random variable.

So I have created an answer here, using NumPy to create the Bernoulli random variable for efficiency (numpy is faster than scipy) and then running the simulation with a stepsize of 0.01 then plotting the solution using matplotlib.

One thing to note that this SDE is one dimensional so we can just store the state and time in separate vectors and plot them at the end.

# Function generating bernoulli trial (your r(t))

def get_bernoulli(p=0.5): 
    '''
    Function using numpy (faster than scipy.stats)
    to generate bernoulli random variable with values 0.5 or 1.6
    '''
    B = np.random.binomial(1, p, 1)
    if B == 0:
        return 0.6
    else:
        return 1.5

This is then used in the simulation as

import numpy as np
import matplotlib.pyplot as plt

dt = 0.01 #step size
x0 = 1# initialize
tfinal = 1 
sqrtdt = np.sqrt(dt)
n = int(tfinal/dt)

# State and time vectors
xtraj = np.zeros(n+1, float)
trange = np.linspace(start=0,stop=tfinal ,num=n+1)

# initialized
xtraj[0] = x0

for i in range(n):
    xtraj[i+1] = xtraj[i] * get_bernoulli(p=0.5)
    
plt.plot(trange,xtraj,label=r'$x(t)$')
plt.xlabel("time")
plt.ylabel(r"$X$")
plt.legend()
plt.show()

Where we assumed the Bernoulli trial is fair, but can be customized to add some more variation.

Simulating an elementary stochastic process in Python

1 Answers1