I need to simulate a huge bunch of compound poisson processes in Matlab on a very fine grid so I am looking to do it most effectively.

I need to do a lot of simulations on the same random numbers but with parameters changing so it is practical to draw the uniforms and normals beforehand even though it means i have to draw a lot more than i will probably need and won't matter much because it will only need to be done once compared to in the order 500*n repl times the actual compound process generation.

My method is the following: Let T be for how long i need to simulate and N the grid points, then my grid is:

t=linspace(1,T,N);

Let nrepl be the number of processes i need then I simulate

P=poissrnd(lambda,nrepl,1); % Number of jumps for each replication
U=(T-1)*rand(10000,nrepl)+1; % Set of uniforms on (1,T) for jump times
N=randn(10000,nrepl); % Set of normals for jump size

Then for replication j:

Poiss=P(j); % Jumps for replication
Uni=U(1:Poiss,j);% Jump times
Norm=mu+sigma*N(1:Poiss,j);% Jump sizes

Then this I guess is where I need your advice, I use this one-liner but it seems very slow:

CPP_norm=sum(bsxfun(@times,bsxfun(@gt,t,Uni),Norm),1);

In the inner for each jump it creates a series of same length as t with 0 until jump and then 1 after, multiplying this will create a grid with zeroes until jump has arrived and then the jump size and finally adding all these will produce the entire jump process on the grid.

How can this be done more effectively?

Thank you very much.