0

Assuming I have a discretized system of SDE of the form

x(:, t+1) = x(:, t) + f1(x(:, t)).*x(:, t)*dt + f2(x(:, t))./(x(:, t).*y(:, t))* sqrt(dt)*rand1;

y(:, t+1) = f2(x(:, t)).*y(:, t)./x(:, t)*dt + f1(x(:, t)).*y(:, t)*sqrt(dt)*rand2;

and I want to simulate the system using 10000 trajectories,

for Time t = 100 days, such that: From Monday to Friday,

f1(x(:, t)) = 2*x(:, t).^2./(y(:, t) + x(:, t) + c), and

f2(x(:, t)) = y(:, t).^2; whereas, Saturdays and Sundays

f1(x(:, t)) = x(:, t)./y(:, t), and f2(x(:, t)) = y(:, t); How can I simulate the SDE system?

Here is my approach

dt = 0.01;
time = 100;
num_iteration = ceil(time / dt);
num_trajectory = 10000; 
%% Initial Values
y0 = 1; 
x0 = 1;
y = zeros(num_trajectory, num_iteration) + y0; 
x = zeros(num_trajectory, num_iteration) + x0; 
days = 0;

for t=1: num_iteration
    current_time = t * dt;
    rand1 = randn(num_trajectory, 1);
    rand2 = randn(num_trajectory, 1);

    if ceil(current_time) == current_time
        days = days+1;

        if (mod(days, 7) | mod(days+1, 7)) == 0
            f1 = 2*x(:, t).^2./(y(:, t) + x(:, t) + c);
            f2 = y(:, t).^2;
        else
            f1 = x(:, t)./y(:, t);
            f2 = y(:, t); 
        end
    end

    x(:, t+1) = x(:, t) + f1*x(:, t)*dt + f2/(x(:, t).*y(:, t))* sqrt(dt)*rand1;
    y(:, t+1) = f2*y(:, t)./x(:, t)*dt + f1*y(:, t)*sqrt(dt)*rand2;   
end
Tasos Papastylianou
  • 21,371
  • 2
  • 28
  • 57
  • 1
    The text does not contain an actual question. What are you trying to achive and what have you done already? – Sven-Eric Krüger May 07 '18 at 12:36
  • For a given `dt`, and assuming that the rows in your variables represent the 10000 "trajectories", then for each row, you require a suitable 'starting point' representing `t=1`, which you will use as `x(:,t)` to find `x(:,t+1)`, i.e. x for t=2. Having t=2 you go on to find x for t=3 etc, up to t = 100. The fact that the `rand` function is involved, introducing some random 'noise' at each step, means that even if the starting point is the same for all rows (e.g. say, x(t=0) = 0 for all rows), the eventual trajectory followed by each row will be different. – Tasos Papastylianou May 07 '18 at 13:52
  • So far this what i have done: – maryam john May 07 '18 at 13:59
  • @TasosPapastylianou, Actually the system is much more complicated than the one i have posted. I did everything you mentioned but i am still not getting the result i am expected to have. – maryam john May 07 '18 at 14:12

1 Answers1

0

Your approach seems fine. You have a logical error in your code though. In the line

if (mod(days, 7) | mod(days+1, 7)) == 0

The expression (mod(days, 7) | mod(days+1, 7)) will always evaluate to 1 (try to figure out why this is), and therefore (mod(days, 7) | mod(days+1, 7)) == 0 will always be false, and your if statement will always pass control to the else part.

Therefore this should instead have been something like

if mod(days, 7) == 0 || mod(days+1, 7) == 0

but this is also confusing (and you haven't documented in code what day of the week day '0' is).

Even clearer would be something like:

if (
  mod (days, 7) == 0 % day is a sunday
  || 
  mod (days, 7) == 6 % day is a saturday
)
    % do stuff
else
    % do other stuff
end 

Better yet, create a small function isWeekend which performs that test for you, leading to super-clear code, like

if isWeekend(days)
  % do stuff
else 
  % do other stuff
end
Tasos Papastylianou
  • 21,371
  • 2
  • 28
  • 57