Построение реализации стохастического процесса на одном графике

Я хочу построить несколько реализаций стохастического процесса в Matlab. Для одной реализации у меня есть следующий код:

N = 80; 
T = dt*N; 
dWt = zeros(1,N);             
S= repmat(S0,1,N);

S(1) = S0;

dWt = sqrt(dt) * randn;

for t=2:N
dWt(t) = sqrt(dt)*randn;
 dSt = k*(mu-S(t-1))*dt + sigma*S(t-1)*dWt(t);
 S(t) = S(t-1)+dSt;

end
 plot(handles.pMeasure, [0:dt:T],[S0,S]);

Я хочу повторить этот цикл n раз и вывести результаты на один график.


matlab plot stochastic-process
person user137425    schedule 07.12.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Вы можете добавить дополнительный цикл for, но лучше всего векторизовать все и вычислить все экземпляры n сразу:

k = ...
mu = ...
sigma = ... 
S0 = ...    % Initial condition
dt = ...    % Time step
n = ...     % Number of instances
N = 80;     % Number of time steps, not counting initial condition
T = dt*N;   % Final time

rng(1);                          % Always seed random number generator
dWt = sigma*sqrt(dt)*randn(n,N); % Calculate Wiener increments
S = zeros(n,N+1);                % Allocate
S(:,1) = S0;                     % Set initial conditions

for t = 2:N+1
    S(:,t) = S(:,t-1) + k*(mu-S(:,t-1))*dt + S(:,t-1).*dWt(:,t-1);
end

plot(handles.pMeasure,0:dt:T,S)

Есть и другие способы оптимизировать это, если хотите, или вы также можете попробовать sde_euler в моем наборе инструментов SDETools Matlab:

k = ...
mu = ...
sigma = ...
dt = ...    % Time step
n = ...     % Number of instances
N = 80;     % Number of time steps, not counting initial condition
T = dt*N;   % Final time

f = @(t,y)k*(mu-y); % Diffusion function
g = @(t,y)sigma*y;  % Drift function
t = 0:dt:T;         % Time vector
S0 = zeros(n,1);    % Initial conditions

opts = sdeset('RandSeed',1,...
              'SDEType','Ito'); % Set random seed, specify Ito SDE
S = sde_euler(f,g,t,S0,opts);   % Simulate

plot(t,S)
person horchler    schedule 07.12.2015
comment
Спасибо! Я получаю сообщение об ошибке. Векторы должны быть одинаковой длины. - person user137425; 08.12.2015
comment
@ user137425: Извините. Это произошло из-за того, что вы определили свой вектор времени с элементами N+1 вместо N. Это должно быть исправлено сейчас. - person horchler; 08.12.2015