psum <- function(vector) {b=vector; b[1]=vector[1]; for (j in 2:length(vector)) b[j]=b[j-1]+vector[j]; b}

gammarw <- function(a,p) 
{unif=runif(10*p,0,1)
 pos=qgamma(unif,a/p,a);
 space=psum(pos);
 time=(1/p)*1:(10*p);
 plot(time,space,
	pch=".",
	sub=paste("Gamma process with shape parameter",a,"and scale parameter",a))}

vgammarw <- function(a,p) 
{unifpos=runif(10*p,0,1)
 unifneg=runif(10*p,0,1)
 pos=qgamma(unifpos,a*a/(2*p),a);
 neg=qgamma(unifneg,a*a/(2*p),a);
 space=psum(pos-neg);
 time=(1/p)*1:(10*p);
 plot(time,space,
	pch=".",
	sub=paste("Variance Gamma process with shape parameter",a*a/2,"and scale parameter",a))}