sim.pp<-function(time,lambda=1){
# Simulates a Poisson process with rate lambda on the
# interval from 0 to time.
sort(time*runif(rpois(1,lambda*time)))
}

sim.cpp<-function(rjump,time,lambda=1){
s<-sim.pp(time,lambda)
y<-rjump(length(s))
list(s=s,y=y)
}


plot.cpp<-function(cpp,time,ylim=range(c(0,cumsum(cpp$y)))){
# Plot a compound Poisson process.
arrivals<-cpp$s
jumps<-cpp$y
x<-c(0,arrivals,time)
y<-c(0,cumsum(jumps),sum(jumps))
plot(x,y,type="s",ylim=ylim,xlab="time",ylab="")
}

jump123<-function(n){sample(3,n,replace=T)}




# Testing

plot.cpp(sim.cpp(jump123,10),10)

plot.cpp(sim.cpp(rexp,10),10)

plot.cpp(sim.cpp(rnorm,10),10)


sim.cpp.rep<-function(rjump,time,nrep,lambda=1){
# Generate a list of nrep Poisson processes.
ans<-list()
for(i in 1:nrep)
ans<-c(ans,list(sim.cpp(rjump,time,lambda)))
ans
}

plot.cpp.rep<-function(cpplist,time){
ylim<-c(0,0)
for(cpp in cpplist)ylim<-range(c(ylim,cumsum(cpp$y)))
for(cpp in cpplist) plot.cpp(cpp,time,ylim=ylim)
}



# Testing

par(mfrow=c(2,2))

clist<-sim.cpp.rep(jump123,10,4)
plot.cpp.rep(clist,10)

clist<-sim.cpp.rep(rexp,10,4)
plot.cpp.rep(clist,10)


clist<-sim.cpp.rep(rnorm,10,4)
plot.cpp.rep(clist,10)


# Now create the plot for class handout.

postscript("sim.compp.ps")
par(mfrow=c(3,2),oma=c(1,1,3,1))

clist<-sim.cpp.rep(jump123,10,6)
plot.cpp.rep(clist,10)

head<-
paste("Compound Poisson Process with Rate 1",
      "and Jumps of 1, 2, or 3",sep="\n")
mtext(head,outer=T,cex=1)


clist<-sim.cpp.rep(rexp,10,6)
plot.cpp.rep(clist,10)

head<-
paste("Compound Poisson Process with Rate 1",
      "and Exponential Jumps",sep="\n")
mtext(head,outer=T,cex=1)


clist<-sim.cpp.rep(rnorm,10,6)
plot.cpp.rep(clist,10)

head<-
paste("Compound Poisson Process with Rate 1",
      "and Normal Jumps",sep="\n")
mtext(head,outer=T,cex=1)

dev.off()


postscript("sim.compp2.ps")
par(mfrow=c(2,2),oma=c(1,1,3,1))

clist<-sim.cpp.rep(jump123,10,4)
plot.cpp.rep(clist,10)

head<-
paste("Compound Poisson Process with Rate 1",
      "and Jumps of 1, 2, or 3",sep="\n")
mtext(head,outer=T,cex=1)


clist<-sim.cpp.rep(rexp,10,4)
plot.cpp.rep(clist,10)

head<-
paste("Compound Poisson Process with Rate 1",
      "and Exponential Jumps",sep="\n")
mtext(head,outer=T,cex=1)


clist<-sim.cpp.rep(rnorm,10,4)
plot.cpp.rep(clist,10)

head<-
paste("Compound Poisson Process with Rate 1",
      "and Normal Jumps",sep="\n")
mtext(head,outer=T,cex=1)

dev.off()