############ Functions ############

# HP Simulation Function
simulate_hawkes <- function(alpha, beta, lambda, n){
  arrivals <- numeric()
  eps <- 1e-6
  # recurrence function
  S <- function(k){
    if(k==1) return(1)
    return( 1 + S(k-1)* exp(-beta * (arrivals[k]- arrivals[k-1])) )
  }
  # the function to be solved to obtain, u
  fu <- function(u, U, k){
    return(log(U) + lambda*(u - arrivals[k]) +
             alpha/beta * S(k) * (1 - exp(-beta * (u - arrivals[k]))))
  }
  fu_prime <- function(u, U , k){
    return( lambda + alpha * S(k) * (exp(-beta * (u - arrivals[k]))))
  }
  # iterative procedure to solve f(u)
  solve_u <- function(U,k) {
    u_prev <- arrivals[k] - log(U)/lambda
    u_next <- u_prev - fu(u_prev, U, k)/ fu_prime(u_prev, U, k)
    while( abs(u_next - u_prev) > eps){
      u_prev <- u_next
      u_next <- u_prev - fu(u_prev, U, k)/fu_prime(u_prev, U, k)
    }
    return(0.5 * (u_prev + u_next))
  }
  #set.seed(1)
  t1 <- -log(runif(1))/lambda
  arrivals <- c(arrivals, t1)
  k <- length(arrivals)
  while(k < n) {
    U <- runif(1)
    t_next <- solve_u( U, k)
    arrivals <- c(arrivals, t_next)
    k <- length(arrivals)
  }
  return(arrivals)
}

# Log Likelihood Function
loglik <- function(params, arrivals){
  alpha_i <- params[1]
  beta_i <- params[2]
  lambda_i <- params[3]
  term_1 <- -lambda_i*arrivals[n]
  term_2 <- sum(alpha_i/beta_i*(exp( -beta_i * (arrivals[n] - arrivals)) - 1))
  Ai <- c(0, sapply(2:n, function(z) {
    sum(exp( -beta_i * (arrivals[z]- arrivals[1:(z - 1)])))
  }))
  term_3 <- sum(log( lambda_i + alpha_i * Ai))
  return(-term_1- term_2 -term_3)
}

# Estimated Intensity function
estimated_intensity <- function(params, arrivals){
  lambda <- params[1]
  alpha <- params[2]
  beta <- params[3]
  Ai <- sapply(1:n, function(z) {
    sum(exp( -beta * (arrivals[z]- arrivals[1:z])))
  })
  return(lambda + alpha *Ai)
}

# Compensator Function
compensator <- function(params, arrivals){
  lambda <- params[1]
  alpha <- params[2]
  beta <- params[3]
  result <- sapply(1:n, function(z){
    lambda*arrivals[z] + sum(alpha/beta*(1-exp(-beta*(arrivals[z]-arrivals[1:z]))))
  })
  return(result)
}
### end of usual Functions ###








# function to find a percentage error as I vary n
errors <- function(n, alpha, beta, lambda){
  # simulate
  sim1 <- simulate_hawkes(alpha, beta, lambda, n)
  # estimate
  MLE <- nlm(loglik, c(2,2,2), hessian = TRUE, arrivals = sim1) 
  # extract estimates
  alphaest<-MLE$estimate[1]
  betaest<-MLE$estimate[2]
  lambdaest<-MLE$estimate[3]
  # percentage error
  alphaerror <- abs(alphaest-alpha)/alpha
  betaerror <- abs(betaest-beta)/beta
  lambdaerror <- abs(lambdaest-lambda)/lambda
  # add errors
  percentageerror <- (alphaerror + betaerror + lambdaerror)/3
  return(percentageerror)
  #return(alphaerror)
  #return(betaerror)
  #return(lambdaerror)
}



# Plot Intensity and Arrivals - Simulation 1
# The only difference between Sim 1 and Sim2 is the "xx = 40 vs 30" in the "ylim=c(4,xx)"
plot_intensity_arrivals1 <- function(sample1, y){
  ## Plot first set of data and draw its axis
  plot(sample1, y,  axes=FALSE,  xlab="", ylab="", type="l",col="black", 
       main="Intensity of the Process vs Arrivals")
  axis(2, ylim=c(0,1),col="black",las=1)  ## las=1 makes horizontal labels
  mtext("Arrivals",side=2,line=2.5)
  box()
  ## Allow a second plot on the same graph
  par(new=TRUE)
  ## Plot the second plot and put axis scale on right
  plot(sample1, estimated_intensity(params, sample1),   
       xlab="", ylab="", ylim=c(4,45), 
       axes=FALSE, type="l", col="red")
  ## a little farther out (line=4) to make room for labels
  mtext("Intensity", side=4, col="red", line=3) 
  axis(4, ylim=c(0,7000), col="red",col.axis="red",las=1) #ylim prob isnt needed here
  ## Draw the time axis
  axis(1)
  mtext("Time",side=1,col="black",line=2.5)  
  ## Add Legend
  legend("topleft",legend=c("Arrivals","Intensity"),
         text.col=c("black","red"),pch=c(16,15),col=c("black","red"))
  # i want to make the black line more visible so i will add it again
  par(new=TRUE)
  plot(sample1, y,  axes=FALSE,  xlab="", ylab="", type="l",col="black")
  #end
}



# Plot Intensity and Arrivals - Simulation 2
# The only difference between Sim 1 and Sim2 is the "xx = 40 vs 30" in the "ylim=c(4,xx)"
plot_intensity_arrivals2 <- function(sample1, y){
  ## Plot first set of data and draw its axis
  plot(sample1, y,  axes=FALSE,  xlab="", ylab="", type="l",col="black", 
       main="Intensity of the Process vs Arrivals")
  axis(2, ylim=c(0,1),col="black",las=1)  ## las=1 makes horizontal labels
  mtext("Arrivals",side=2,line=2.5)
  box()
  ## Allow a second plot on the same graph
  par(new=TRUE)
  ## Plot the second plot and put axis scale on right
  plot(sample1, estimated_intensity(params, sample1),   
       xlab="", ylab="", ylim=c(4,40), 
       axes=FALSE, type="l", col="red")
  ## a little farther out (line=4) to make room for labels
  mtext("Intensity",side=4,col="red",line=3) 
  axis(4, ylim=c(0,7000), col="red",col.axis="red",las=1) #ylim prob isnt needed here
  ## Draw the time axis
  axis(1)
  mtext("Time",side=1,col="black",line=2.5)  
  ## Add Legend
  legend("topleft",legend=c("Arrivals","Intensity"),
         text.col=c("black","red"),pch=c(16,15),col=c("black","red"))
  # i want to make the black line more visible so i will add it again
  par(new=TRUE)
  plot(sample1, y,  axes=FALSE,  xlab="", ylab="", type="l",col="black")
}




# multiple HP simulations then return the average of the 10 sims
multsim <- function(alpha, beta, lambda, n){
  sim_1 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_2 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_3 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_4 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_5 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_6 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_7 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_8 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_9 <- simulate_hawkes(alpha, beta, lambda, n)
  sim_10<- simulate_hawkes(alpha, beta, lambda, n)
  #average sim
  avgsim = c()
  for (i in 1:n){
    avgsim[i] = (sim_1[i]+sim_2[i]+sim_3[i]+sim_4[i]+sim_5[i]+sim_6[i]+sim_7[i]+sim_8[i]+sim_9[i]+sim_10[i])/10
  }
  return(avgsim)
}





#end - STOP