2

我目前有模拟几何布朗运动的代码，由http://www-math.bgsu.edu/~zirbel/sde/matlab/index.html提供。

但是，我想生成 1,000 个模拟并将它们显示在图表中。

我目前生成单个模拟的代码如下：

% geometric_brownian(N,r,alpha,T) simulates a geometric Brownian motion 
% on [0,T] using N normally distributed steps and parameters r and alpha

function [X] = geometric_brownian(N,r,alpha,T)

t = (0:1:N)'/N;                   % t is the column vector [0 1/N 2/N ... 1]

W = [0; cumsum(randn(N,1))]/sqrt(N); % S is running sum of N(0,1/N) variables

t = t*T;
W = W*sqrt(T);

Y = (r-(alpha^2)/2)*t + alpha * W;

X = exp(Y);

plot(t,X);          % plot the path
hold on
plot(t,exp(r*t),':');
axis([0 T 0 max(1,exp((r-(alpha^2)/2)*T+2*alpha))])
title([int2str(N) '-step geometric Brownian motion and its mean'])
xlabel(['r = ' num2str(r) ' and alpha = ' num2str(alpha)])
hold off
4

2 回答 2

6

That code cannot be used directly to simulate 1,000 paths/simulations. Unfortunately, it has not been vectorized. The easiest way to do what you want is to use a for loop:

N = 1e3;
r = 1;
alpha = 0.1;
T = 1;
npaths = 1e3;          % Number of simulations

rng(0);                % Always set a seed
X = zeros(N+1,npaths); % Preallocate memory
for i = 1:n
    X(:,i) = geometric_brownian(N,r,alpha,T);
    hold on
end
t = T*(0:1:N).'/N;
plot(t,exp(r*t),'r--');

This is rather slow and inefficient. You will need to modify the function a lot to vectorize it. One thing that would improve performance is if you at least removed the plotting code from inside the function and ran that separately after the loop.

Another alternative might be to use the sde_gbm function in my SDETools toolbox, which is fully-vectorized and much faster:

N = 1e3;
r = 1;
alpha = 0.1;
T = 1;
npaths = 1e3;        % Number of simulations

t = T*(0:1:N)/N;     % Time vector
y0 = ones(npaths,1); % Vector of initial conditions, must match number of paths
opts = sdeset('RandSeed',0,'SDEType','Ito'); % Set seed
y = sde_gbm(r,alpha,t,y0,opts);

figure;
plot(t,y,'b',t,y0*exp(r*t),'r--');
xlabel('t');
ylabel('y(t)');
title(['Geometric Brownian motion and it's mean: ' int2str(npaths) ...
       ' paths, r = ' num2str(r) ', \alpha = ' num2str(alpha)]);

In either case, one obtains a plot that looks something like this

               enter image description here

于 2013-09-08T20:38:42.340 回答
0

要执行1000 次模拟，直接的方法是：

Nsims = 1000;
N=10^15;         % set to length of individual sim                    
r = 1;
alpha = 0.1;
T = 1;

t = (0:1:N)'/N;                   
t = (T*(r-(alpha^2)/2))*t;
W = cat(1,zeros(1,Nsims),cumsum(randn(N,Nsims))); 
W = W*(sqrt(T)*alpha/sqrt(N));
Y = repmat(t,1,Nsims) + W;
X = exp(Y);

绘图就像以前一样

plot(t,X);             % plots ALL 1000 paths
%   plot(t,X(:,paths));   % use instead to show only selected paths (e.g. paths =[1 2 3])
hold on
plot(t,exp(r*t),':');
axis([0 T 0 max(1,exp((r-(alpha^2)/2)*T+2*alpha))])
title([int2str(N) '-step geometric Brownian motion and its mean'])
xlabel(['r = ' num2str(r) ' and alpha = ' num2str(alpha)])
hold off

对于相对较短（小）的模拟集，循环遍历您的代码或执行上述操作应该可以。对于重型模拟，您可能会受益于 Horchler 承诺的速度优势。

于 2013-09-08T15:02:40.250 回答