Это больше похоже на формулу частичной интеграции, или эквивалентного правилу продукта
(f*b)'=f'*b+f*b' or
d(f*b)=df*b+f*db = f'*b*dt+f*db
дифференциации, чем любой стохастического дифференциального уравнения. В этом смысле там ничего не решить, за исключением случаев, когда вы хотите количественно проверить правильность правила с b
, являющимся винеровским процессом.
Или вы хотите проверить ошибку численных методов. Но тогда первым шагом было бы определить численный метод, который, по-видимому, является методом Эйлера или Эйлера-Маруямы.
См wikipedia: Euler-Maruyama method, wikipedia: Milstein method, и очень читаемый P. Forsyth: An introduction to Computational Finance without Agonizing Pain.
Для процесса Wiener b
вам нужно будет держать массив приращений db=randn(1,N)*sqrt(dt)
(и сохранить корень из f=sin(t)
). Тогда левая часть представляет собой сумму по f.*db
и на правой стороне вы должны были бы массив частичных сумм по db
, которая в теории должна быть вычислена как
b(i+1)=b(i)+db(i)
Использование более целесообразный функции MatLab
b = cumsum(db)
предоставит вам смещенную версию, где b(i)=b(i-1)+db(i)
. В общей сумме изображения этот сдвиг может привести к ошибке величиной sqrt(dt)
, что было бы важно, если вы хотите ограничить ошибку величиной dt
.
После этой операции действуют следующие условия: f(N)*b(N)-f(0)*b(0)
и сумма свыше cos(t(i))*b(t(i))*dt(i)
, i=0,...,N-1
.
Таким образом, вы сравните с левой
cumsum(sin(t).*db)
и справа
sin(t).*b - cumsum(cos(t).*b).*dt
и это должно дать практически одинаковые значения.
В кузена MatLab в Scilab следующий сценарий дает совершенно перекрывающихся результаты:
N=6000; dt=2*%pi/N;
t=0:1:N; t=t*dt;
dW=rand(t,"normal")*sqrt(dt);
W = cumsum(dW);
f=sin(t); fp=cos(t);
ex1=cumsum(f.*dW);
ex2=f.*W;
ex3=cumsum(fp.*W).*dt;
clf;
subplot(211);
plot(t+5*dt,ex1,t,ex2-ex3)
subplot(212);
plot(t,ex2-ex1-ex3)
можно увидеть идентичность также непосредственно в дискретизации пути применения индекса сдвигов. Установка b=cumsum(db)
, f=sin(t)
и fp(t)=cos(t)
производной элемент n
из cumsum(f.*db)
является и преобразуется как:
sum(i=1..n) f(i)*db(i) = sum(i=1..n) f(i)*(b(i)-b(i-1))
= -f(1)*b(0)+sum(i=1..n) (f(i)-f(i+1))*b(i) + f(n+1)*b(n)
= f(n)*b(n) - sum(i=1..n) fp(i)*b(i)*dt
+ (f(n+1)-f(n))*b(n) - sum(i=1..n) (f(i+1)-f(i)-fp(i)*dt)*b(i)
с виртуальным значением b(0)=0
. Следующая последняя строка содержит n
-й элемент f.*b-cumsum(fp.*b)*dt
, а в последней строке - условия ошибки дискретизации с величинами dt
и (t(n)-t(1))*dt
.
StackOverflow не поддерживает LaTeX, поэтому я взял ваш синтаксис, отобразил их как изображения и встроил их в сообщение. – rayryeng
fprime (s) = cos (s), b (s) = sum (randn (N * s/t, 1) * sqrt (dt)) и ds = dt ... вы должны использовать одни и те же случайные числа для правая и левая стороны, естественно –