How to compute the mean of several stochastic realizations

I have a code that can be used to plot sample paths of the stochastic differential equation competition model

clear
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;
for j=1:num_realizations
    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;
    end
    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

My question is how to plot the average of each population, including the initial conditions?

[John D'Errico]

Simulate it for a long enough time, and use mean. Where is the problem? Are you asking to compute an analytical mean in some way?

[Fares]

Hi John, thanks for replying!

I know how to calculate means for simple examples but I am not sure how to do it when there is a loop like the case above.

[Fares]

I am basically looking for the mean of x1(for example) as an array where the first element represents the average of the first elements of all realizations, second element represents the average of the second elements of all realizations and so on.

[Umar]

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.

[Fares]

Hi Umar, thank you so much! This is very helpful.

The two arrays x1 and x1_mean has two different dimensions. One is a row vector and the other is a column vector. This is easy to fix. We can change from zeros(k,1) into zeros(1,k). However, the above code is missing one value becuase array x1 is of length 1X5001 but the mean is of length 1x5000. Shoud we change from zeros(1,k) to zeros(1,k+1) to account for the missing value?

[Umar]

Hi Fares,

Thankyou for noticing that. I do apologize, changing zeros(k,1) to zeros(1,k) is a step in the right direction. However, you highlighted that there is a missing value due to the difference in lengths between x1 (1x5001) and x1_mean (1x5000). In order to account for this missing value, adjusting from zeros(1,k) to zeros(1,k+1) could indeed be a viable solution.

[Fares]

Thank Umar for your response. In case I want also to compute the variance, do I need to re-run simulations?

[Umar]

Hi Fares,

To answer your question and save time and computational resources, I would suggest computing the variance without re-running simulations. To accomplish this task, you can utilize the mean values already calculated in your existing simulation. After calculating the mean values for each time step, you can then compute the variance using the formula for variance, let me illustrate by example,

% Calculate Variance of x1

x1_variance = var(x1_mean);

disp('Variance of x1:');

disp(x1_variance);

So, using using the var() function in Matlab on the array x1_mean that stores the mean values, you can efficiently compute the variance without the need to re-run the simulations. You probably familiar with var() function, in case you seek guidance, please refer to

https://www.mathworks.com/help/matlab/ref/var.html?searchHighlight=var&s_tid=srchtitle_support_results_1_var

Please let me know if you have any further questions.

[Fares]

Thanks Umar for your valuable support. I am not looking for the variance of the mean of x1. I am looking for the variance of x1 as an array again where the first element represents the variance of the first elements of all realizations, second element represents the variance of the second elements of all realizations and so on. My goal is to to compute the standard deviation so I can plot all trajectories within +/- one standartds deviation.

[Umar]

Hi Fares,

Thanks for clarifying this. You can define the realizations array containing multiple realizations. Then, ise the var function with your defined parameters to calculate the variance along columns and store it in variance_array. Afterwards, compute the standard_deviation by taking the square root of the variance array. As far ploting Trajectories, you need to calculate the upper_bound and lower_bound arrays by adding and subtracting the standard deviation from the realizations array, respectively. Then,create a plot to visualize the trajectories within +/- one standard deviation.

Hope these steps will help you to achieve your goal now.

[Fares]

Thanks Umar for your help. Really appreciated!

[Umar]

No problem Fares, glad to help out.

Answer 1

[Umar]

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.