# Renewal process simulation 
# STAT 753
# April 14, 2020


##############################
# Simulate a renewal process 
##############################

############################################################
# Version 1 uses a uniform distribution of interevent times 
# Default is U(0,1), i.e. min=0 and max=1
############################################################

library(markovchain)

renewal.uniform.sim <- function (min=0, max=0, num.events, num.sims = 1, t0 = 0) 
{
  if (num.sims == 1) {
    x = t0 + cumsum(runif(n = num.events, min, max))
    return(c(t0,x))
  }
  else {
    xtemp = matrix(runif(n = num.events * num.sims, min, max), num.events)
    x = t0 + apply(xtemp, 2, cumsum)
    return(rbind(rep(t0, num.sims), x))
  }
}


# Run the function above, num.sims default is 1
renewal.uniform.sim(min=0, max=1, num.events = 10)
renewal.uniform.sim(min=0, max=1, num.events = 10, num.sims = 5)


# Run the function and plot the time series (for 1 simulation)
num.events = 100
run1 <- renewal.uniform.sim(0, 1, num.events)
plot(run1, 0:100, xlab = "Time", ylab = "Number of Events", 
     main = "Renewal Process Sample Path: Uniform Interarrival Distribution", pch = 19, col="blue")

head(run1)
tail(run1)

# Compare to a Poisson process with the same mean interarrival time
## Exercise 


#########################################################
# Version 2 uses a beta distribution of interevent times 
# Default is beta(a=2,b=2) -- normal shaped 
#########################################################

renewal.beta.sim <- function (a=2, b=2, num.events, num.sims = 1, t0 = 0) 
{
  if (num.sims == 1) {
    x = t0 + cumsum(rbeta(n = num.events, a, b))
    return(c(t0,x))
  }
  else {
    xtemp = matrix(rbeta(n = num.events * num.sims, a, b), num.events)
    x = t0 + apply(xtemp, 2, cumsum)
    return(rbind(rep(t0, num.sims), x))
  }
}


# Run the function above, num.sims default is 1, uses default a and b values
renewal.beta.sim(num.events = 10)

# Change the values of the beta parameters
# Interesting cases: a=b=0.5 -- U shaped; a=2, b=5 -- skewed left; a=1, b=3 -- decreasing exponential-like
a = 0.5
b = 0.5
renewal.beta.sim(a, b, num.events = 10)
run <- renewal.beta.sim(a, b, num.events = 10, num.sims = 5)

run2 <- renewal.beta.sim(a, b, num.events = 100)
plot(run2, 0:100, xlab = "Time", ylab = "Number of Events", 
     main = "Renewal Process Sample Path: Beta Interarrival Distribution", pch = 19, col="tomato")


# Run the function and plot the time series (for 1 simulation)
a = 1
b = 3
num.events = 100
run3 <- renewal.beta.sim(a, b, num.events)
plot(run3, 0:100, xlab = "Time", ylab = "Number of Events", 
     main = "Renewal Process Sample Path: Beta Interarrival Distribution", pch = 19, col="darkgreen")