How to compute the mean of several stochastic realizations

Hi Fares,

Here is the modified code snippet with the mean calculation implemented:

>> a10 = 2; a20 = 1.5; a11 = 0.03; a12 = 0.02; a21 = 0.01; a22 = 0.04; x1(1) = 50; x2(1) = 25; k = 5000; T = 5; dt = T / k; num_realizations = 10; x1_mean = zeros(k, 1); % Initialize array to store mean of x1

for j = 1:num_realizations x1_realization = zeros(k, 1); % Initialize array to store x1 values for each realization for i = 1:k rn = randn(2, 1); f1 = x1(i) * (a10 - a11 * x1(i) - a12 * x2(i)); f2 = x2(i) * (a20 - a21 * x1(i) - a22 * x2(i)); g1 = x1(i) * (a10 + a11 * x1(i) + a12 * x2(i)); g2 = x2(i) * (a20 + a21 * x1(i) + a22 * x2(i)); x1(i + 1) = x1(i) + f1 * dt + sqrt(g1 * dt) * rn(1); x2(i + 1) = x2(i) + f2 * dt + sqrt(g2 * dt) * rn(2); x1p = [x1(i + 1) > 0]; x2p = [x2(i + 1) > 0]; x1(i + 1) = x1(i + 1) * x1p; x2(i + 1) = x2(i + 1) * x2p;

        x1_realization(i) = x1(i); % Store x1 values for each realization
    end
    x1_mean = x1_mean + x1_realization; % Accumulate x1 values for mean calculation

    plot([0:dt:T], x1, 'Color', rand(1, 3), 'LineWidth', 1);
    hold on
    plot([0:dt:T], x2, 'Color', rand(1, 3), 'LineWidth', 1);
    ylabel('Population Size'); xlabel('Time');
end

x1_mean = x1_mean / num_realizations; % Calculate mean of x1 disp('Mean of x1 for each time step:'); disp(x1_mean);

x1_mean array is used to accumulate the values of x1 from each realization. At the end of all realizations, the mean of x1 at each time step is calculated by dividing the accumulated values by the total number of realizations. This implementation will provide you with an array x1_mean where each element represents the average of x1 values across different realizations at that specific time step. This should get you started with rest such as second element representing the average of the second elements of all realizations and so on. Feel free to customize the code based on your needs.

Please see attached results.