function [Sp Sq] = SimCRR(T, N, mu, sigma, r);

% tim step size
dt = T/N;

% branching probabilities
p = 0.5* ( 1 + (mu-0.5*sigma^2)/sigma * sqrt(dt) );
q = 0.5* ( 1 + (r-0.5*sigma^2)/sigma * sqrt(dt) );

% generate Bernoullis
u = rand(1,N);
xp = 2*(u < p) - 1;
xq = 2*(u < q) - 1;

Sp(1) = 1;
Sp(2:N+1) = Sp(1) * exp( sigma*sqrt(dt) * cumsum(xp) );

Sq(1) = 1;
Sq(2:N+1) = Sq(1) * exp( sigma*sqrt(dt) * cumsum(xq) );



end