function [a1 a2]=w(M,N) % solving dX = alpha(a-X)dt + sqrt(X)sigma dW % with a=5 alpha=1 sigma=0.3 X(0)=1 % to find E[X(2)] % using M simulations and N subdivisions of [0,2] a=5; alpha=1; sigma=0.3; X0=1; T=2; h=T/N; r=randn(M,N); % normal random numbers X=X0*ones(M,1); for i=1:N X = X+h*alpha*(a-X)+sqrt(h)*sigma*sqrt(X).*r(:,i); end a1=mean(X); a2=std(X);
Running
[a b]=w(10000,20)gave answers a=4.5221 and b=0.4337. The stochastic error is plus or minus 0.004. Further runs gave a=4.5184, 4.5159, 4.5147, 4.5060, 4.5113, 4.5168.
To implement for different values of h between 0.1 and 0.3 I did the following:
h=[0.1:0.01:0.3]';
hs=zeros(size(h));
as=zeros(size(h));
bs=zeros(size(h));
for i=1:size(h,1)
N=floor(2/h(i));
[as(i) bs(i)]=w(100000,N);
hs(i)=2/N
end
The b's varied little: from 0.4383 to 0.4744, giving a stochastic error
of about plus or minus 0.0015. Doing plot(hs,as)and doing a linear fit gave
Clearly the error behaves linearly with h, as expected (the stochastic error is relatively small compared with the error due to the use of large values of h).
Back to main course page
Back to my main teaching page
Back to my home page