/* Gamma Process */
new; cls;
alpha=10;
beta=2;
X0=50;
T=1;
n=252;
X=X0|gamp(alpha,beta,X0,T,n);

library pgraph;
graphset;
xy(seqa(0,T/n,n+1),X);




proc gamp(alpha,beta,X0,T,n);
   local delt,G,X,i;
   delt=T/n;
   G=rndgamma(alpha*delt,beta,n);
   X=zeros(n,1);
   X[1]=X0+G[1];
   i=1;
   do while i<=n-1;
      X[i+1]=X[i]+G[i+1];
      i=i+1;
   endo;
   retp(X);
endp;

proc rndgamma(alpha,beta,nn);
   local x;
   if alpha<=0 or beta<=0;
      errorlog "ERROR: Parameter alpha and beta must be positive.";
      retp(".");
   endif;
   x=rndu(nn,1);
/* U(0,1) to Gamma(alpha,beta) */
   x=beta*gammaii(alpha*ones(nn,1),x);
   retp(x);
endp;