new;

library pgraph;

proc (3) = saut(x0,mu,sigma,kappa,lambda,t0,TT,N,Ns);
  local h,t,xt,u,dN,i,xi,ti,dN_,lambda_;
  local mu:proc,sigma:proc,kappa:proc,lambda: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);
  dN = zeros(N+1,Ns);

  i = 1;
  do until i > N;
    xi = xt[i,.]; ti = t[i];
    lambda_ = lambda(ti,xi);
    dN_ = rndu(1,Ns) .< (lambda_*h);
    xt[1+i,.] = xi + mu(ti,xi)*h
                   + sigma(ti,xi).*u[i,.] + kappa(ti,xi).*dN_;
    dN[1+i,.] = dN_;
    i = i + 1;
  endo;
  retp(t,xt,dN);
endp;



x0 = 2;

proc mu(t,x);  		retp(0.0002*x);endp;
proc sigma(t,x);  	retp(0);	endp;
proc kappa(t,x);  	if t <= 2;   	retp(.01); 
	else;	   	retp(5.);  endif;	endp;

rndseed 123;
proc lambda(t,x);   
	if t <= 2;   	retp(5); 
	else;	   	retp(0.5);  endif;	endp;
{t,x2,dN} = saut(x0,&mu,&sigma,&kappa,&lambda,0,10,100,1);

xy(t,x2);