# Simulation - Code



##### Simulation #####

# Simulate n = 500
n <- 500
alpha <- 4
beta <- 5
lambda <- 0.5
params = c(lambda, alpha, beta) #needed for some functions, not simulation function
# now call the function (alpha, beta, lambda, n)
# name the vector of arrivals as sim1
sim1 <- simulate_hawkes(alpha, beta, lambda, n)
# plot the simulated counting process
y<-c(1:n)
plot(sim1, y,  type="l",  main="Hawkes Process Simulation (n = 500)", ylab="Arrivals", xlab="Time")
# MLE using LogLike (Initial guesses in blue)
MLE <- nlm(loglik, c(2,2,2), hessian = TRUE, arrivals = sim1) 
MLE
# print MLE
paste( c("alpha", "beta", "lambda"), round(MLE$estimate,2), sep=" = ")
# guess:  "alpha = 3.54"  "beta = 4.36"   "lambda = 0.57"
# answer: "alpha = 4.00"  "beta = 5.00"   "lambda = 0.50"
alphaest<-MLE$estimate[1]
betaest<-MLE$estimate[2]
lambdaest<-MLE$estimate[3]
# Branching Ratio
Branch=alphaest/betaest
Branch
#
# simulate a HP using est. params
sim2 <- simulate_hawkes(alphaest, betaest, lambdaest, n) #rename
plot(sim2, y,  type="l",  main="Hawkes Process Simulation (based on estimated parameters)", 
     ylab="Arrivals", xlab="Time")
#plot both together
plot(sim1, y,  type="l",  main="Hawkes Process Simulation vs Resimulated Process", 
     ylab="Arrivals", xlab="Time")
lines(sim2, y, col="blue")
legend("topleft",legend=c("Original Simulation","Resimulation from Estimated Parameters"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))
# we can see its the same distribution but at a bit of a lag
plot_intensity_arrivals1(sim1, y)
# plot residuals
diff<-sim2-sim1
qqnorm(diff, pch = 1, frame = FALSE)
qqline(diff, col = "steelblue", lwd = 2)


# Simulate n = 1000 (plus new parameters)
n <- 1000
alpha <- .225
beta <- 3
lambda <- 1
params = c(lambda, alpha, beta)
# now call the function (lambda, alpha, beta, n)
sim1<-simulate_hawkes(alpha, beta, lambda, n)
# plot the simulated counting process
y<-c(1:n)
plot(sim1, y,  type="l",  main="Hawkes Process Simulation (n = 1000)", ylab="Arrivals", xlab="Time")
# MLE using LogLike (Initial guesses in blue)
MLE <- nlm(loglik, c(1,1,1), hessian = TRUE, arrivals = sim1) 
MLE
# print MLE
paste( c("alpha", "beta", "lambda"), round(MLE$estimate,2), sep=" = ")
# guess:  "alpha = 2.90"  "beta = 3.87"   "lambda = 1.06"
# answer: "alpha = 3.00"  "beta = 4.00"   "lambda = 1.00"
alphaest<-MLE$estimate[1]
betaest<-MLE$estimate[2]
lambdaest<-MLE$estimate[3]
# Branching Ratio
Branch=alphaest/betaest
Branch
#
# simulate a HP using est. params
sim2 <- simulate_hawkes(alphaest, betaest, lambdaest, n) #rename
plot(sim2, y,  type="l",  main="Hawkes Process Simulation (based on estimated parameters)", 
     ylab="Arrivals", xlab="Time")
#plot both together
plot(sim1, y,  type="l",  main="Hawkes Process Simulation vs Resimulated Process", 
     ylab="Arrivals", xlab="Time")
lines(sim2, y, col="blue")
legend("topleft",legend=c("Original Simulation","Resimulation from Estimated Parameters"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))
# we can see by doubling the sample size it reduces the drift
plot_intensity_arrivals2(sim1, y)
# plot residuals
diff<-sim2-sim1
qqnorm(diff, pch = 1, frame = FALSE)
qqline(diff, col = "steelblue", lwd = 2)