new;
library pgraph;

proc (2) = mbg(x0,mu,sigma,t0,TT,N,Ns);
  local h,t,xt,u,i,xi,ti;
  local mu:proc,sigma:proc;

  h = (TT-t0)/N;
  t = seqa(t0,h,N+1);

  xt = zeros(N+1,Ns);
  xt[1,.] = x0.*ones(1,Ns);
  u = rndn(N,Ns)*sqrt(h);

  i = 1;
  do until i > N;
    xi = xt[i,.]; ti = t[i];
    xt[1+i,.] = xi + mu(ti,xi)*h + sigma(ti,xi).*u[i,.];
    i = i + 1;
  endo;
  retp(t,xt);
endp;

x0 = 25;
rndseed 123;	
proc mu(t,x);  		retp(0.0002*x);endp;
proc sigma(t,x);  	retp(.09*x);	endp;
{t,x1} = mbg(x0,&mu,&sigma,0,10,100,1000);

/*graphset;
  fonts("simplex simgrma");
  _pdate = ""; _pnum = 2;
  title("Mouvement Brownien Geometrique : "\
        "dX]t[ = \202m\201X]t[dt + \202s\201X]t[dW]t["\)
*/
xy(t,x1);