simulate N trajectories for each of two independent geometric brownian motion

Hello,

I have to wirte a code in matlab to simulate N trajectories of steel spot prices and N trajectories of electricity spot prices which are governed by independent geometric Brownian motions. i have already wrote a code to generate N trajectories for the first stock, but i don't know how to do it for both of them. can anyone help me plss..??

thanks in advance

my code is:

 %Simulate geometric Brownian motion in [t0,t0+T]
    %INPUT:
       % nb_traj  : Number of trajectories
       % N  : Number of time steps
       % t0 : initial time
       % S1Begin : initial value 
       % T  : Length of the simulation interval [t0,t0+T]
       % drift , sigma : parameters of the process
    %OUTPUT:
        %t  : Times at which the trajectory is monitored
        %S  : Simulation of nb_traj trajectories of geometric Brownian motion

T = 10; %time to maturity S1Begin = 200; %Initial price of steel sigma = 0.627; %Volatility of steel drift = 0.025; %Drift, risk-neural of steel nb_traj = 100; %number of trajectories dt = 1/25; %time steps S1Optimal = 30; %the trigger for undertaking the investment

[S, t] = geom_brownian( sigma, T, S1Begin, drift, dt, nb_traj);

steps = T/dt + 1; tInf = zeros(1,nb_traj); for j = 1:nb_traj; tInf(j) = T; end for tr = 1:nb_traj; for i = 1:steps; if ((S(i,tr) <= S1Optimal) && t(i) < (tInf(tr))) tInf(tr) = t(i); end end end

sumTimes = 0; for k = 1:nb_traj; sumTimes = sumTimes + tInf(k); end timeToInvest = sumTimes/nb_traj;

figure('Color',[0.9412 0.9412 0.9412 ]); plot(t,S, 'linewidth',2); axis([t(1) t(end) min(S(:)) max(S(:))]); xlabel('time');ylabel('value S'); title(sprintf('%d trajectories', size(S,2)));

function [S, t] = geom_brownian( sigma, T, S1Begin, drift, dt, nb_traj) % Geomteric Brownian Motion % uso: % S = geom_brownian( sigma, T, S0, mu, step) % esempio: % [S, t] = geom_brownian( .12, 5, 50, 100, .05, 1/255);

nT = ceil(T/dt); % numero di passi = T/step W = sigma * sqrt(dt) * cumsum(randn(nT,nb_traj)); % W_k = sigma*sqrt(step)*(U_k + ... + U_0) c = repmat((drift - sigma^2/2) *dt * (1:nT)',1,nb_traj); S = [repmat(S1Begin,1,nb_traj); S1Begin * exp( c + W)]; %S = [repmat(S0,1,nb_traj); S0 * exp( c + W)]; if nargout > 1 t = [0;dt * (1:nT)'];

end