Я хочу смоделировать траекторию наверняка стохастического процесса в непрерывном времени, используя Python. Я написал следующую функцию, чтобы имитировать это (конечно, np
относится к numpy
):
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)
Функция independent_brownian
просто создает путь для стандартного броуновского движения. Для полноты вот это:
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)
Бывает, что математика, стоящая за моей моделью, подразумевает, что процесс $V_t$, представленный его дискретизацией V
в приведенном выше коде, должен быть положительным. Однако это может быть случай, когда массив dB
содержит отрицательные элементы, большие по абсолютной величине. Я хотел бы автоматизировать следующую процедуру «выбора»:
simulate_V
;V[i]
упадет ниже нуля, прервите процесс, попробуйте другую последовательность dB
и начните сначала;V
положительны;Каков хороший способ автоматизировать эту процедуру? Прямо сейчас я понимаю, что если V[i]
опустится ниже нуля, я получу nan
в numpy, но он не выдаст ошибку или не остановит процесс.
Ради общности и поскольку я не очень хорошо знаком с броуновским движением, я реализую общий механизм для ситуации, в которой у нас есть переменная dB
, которую мы используем для создания другой переменной V
, и нам нужно повторить эту генерацию из V
, пока не будет выполнено определенное условие.
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)
Это дало мне следующий вывод:
[2.27582634 0.77008881 0.28388536]
[2.55101104 3.10944337 0.55829105]
И вы можете видеть, что условие на V
выполняется.
Если необходимы какие-либо другие разъяснения, дайте мне знать.