2015-03-26 5 views
1

Я пытаюсь приблизить решение:Аппроксимируя Решение стохастического интеграла

где и для обеих сторон этого уравнения.

Мой код для левой стороны:

N = 2000; Tend = 2*pi; dt = Tend/N; t = 0:dt:Tend; 
f = sin(t)*sqrt(dt); 
f = [0 ff(1:end-1)]; 
[fL,junk] = meshgrid(f,1); 
dW = cumsum([0 randn(1,N)].*fL,2); 

, но я не могу понять, с правой стороны на все, что гораздо сложнее. Может ли кто-нибудь помочь?

+1

StackOverflow не поддерживает LaTeX, поэтому я взял ваш синтаксис, отобразил их как изображения и встроил их в сообщение. – rayryeng

+0

fprime (s) = cos (s), b (s) = sum (randn (N * s/t, 1) * sqrt (dt)) и ds = dt ... вы должны использовать одни и те же случайные числа для правая и левая стороны, естественно –

ответ

0

Это больше похоже на формулу частичной интеграции, или эквивалентного правилу продукта

(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.

+0

Хм ... может быть, я должен спросить, в чем разница, когда я генерирую db и просто b, где b - это процесс wiener. Эта проблема дала понять, что я не думаю, что больше понимаю эту разницу. – Islands

+0

db - белый шум, b - процесс Винера, также известный как броуновское движение (физическое броуновское движение - процесс Орнштейна-Уленбека, есть некоторые опасения в схеме именования). Если 'dt (i) = t (i + 1) -t (i)' then 'db (i) = randn() * sqrt (dt (i))' и 'b (i + 1) = b (i) + дБ (я) '. – LutzL

+0

Итак, b (0) = 0? – Islands