# PCS - Code

# Predicted Cascade Size = PCS
PCS <- function(lambda,beta,alpha,T,tao,omega){
  # create required variables
  xi=alpha/beta
  r=omega-T
  # create the "a" vector
  a=c()
  n=length(tao)
  for (i in 1:n){
    a[i]=T-tao[i]
  }
  sum=0
  m=0
  # loop to find the second part of the function
  for (i in 1:(n-1)){
    newsum=(xi*exp(-beta*a[i]))/(1-xi)*(1-exp(-beta*(1-xi)*r))
    sum=sum+newsum                                               
  }
  # add the second part of the function to the rest
  m = n + (lambda/(1-xi))*(r+xi/(beta*(1-xi))*(exp(-beta*(1-xi)*r)-1)) + sum
  return(m)
}
# top of page 13 explains a lot of the letters
# this function is made to calculate eq 54 page 16
# https://arxiv.org/pdf/2001.09490.pdf








#####################################################
#####################################################
############## Testing Joey's Function ##############
#####################################################
#####################################################











############################################################################################
### ##############################        n = 500        ############################### ###
############################################################################################

# test joeys function on n=500 simulation of HP
n <- 500
alpha <- 4
beta <- 5
lambda <- 0.5
params = c(lambda, alpha, beta) #needed for some functions, not simulation function
# 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")
sim1

# alpha, beta, lambda all defined above
# I show the function 100 out of the 500 arrivals
T=sim1[100]       # T is the time of the last observed arrival 
tao=sim1[1:100]   # tao is the vector of arrival times observed
omega=sim1[500]   # omega is the time at which i will predict how many arrivals will have come in total

PCSend<-PCS(lambda,beta,alpha,T,tao,omega)
PCSend
# expecting m = 500 as a predicted cascade size
# got 472.4123
# very close
# here, i only let the function see the first 100 arrivals, and it predicted the next 400!!!!!!

#plot this
plot(sim1, y,  type="l",  main="Hawkes Process Simulation with Predicted Cascade Size (n = 500)", 
     ylab="Arrivals", xlab="Time")
PCSline = c(100,PCSend)
PCSy=c(T,omega)
lines(PCSy, PCSline, col="blue")
legend("topleft",legend=c("Original Simulation","Predicted Cascade Size (Estimated after arrival 100)"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))









############################################################################################
### ##############################       n = 1000        ############################### ###
############################################################################################


#testing joeys function with n=1000
n <- 1000
alpha <- 5
beta <- 7
lambda <- 1.5
params = c(lambda, alpha, beta) #needed for some functions, not simulation function
# 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")
sim1

# alpha, beta, lambda already defined above
T=sim1[500]
tao = sim1[1:500]
omega=sim1[1000]
#omega=sim1[1000] #this is the same as omega=182....

PCS(lambda,beta,alpha,T,tao,omega)
#expecting m = 1000 as a predicted cascade size
# got 1013.94
# extremely close
# only showed the function half of what it was working with

# plot this - this way of plotting is way more complicated than it needs to be
time=seq(500, 1000, by=10)
mvector=c()
for (i in time){
  omega=sim1[i]
  m=PCS(lambda,beta,alpha,T,tao,omega)
  m
  mvector=c(mvector,m)
}
mvector

#make y axis
y=c()
for (i in time){
  ynew=sim1[i]
  y=c(y,ynew)
}
y

#plot the predicted over the actual
yold=c(1:n)
plot(sim1, yold,  type="l",  main="Hawkes Process Simulation with Predicted Cascade Size (n = 1000)", 
     ylab="Arrivals", xlab="Time")
lines(y, mvector , col="blue")
legend("topleft",legend=c("Original Simulation","Predicted Cascade Size (Estimated after arrival 500)"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))









############################################################################################
### ##############################       n = 2000        ############################### ###
############################################################################################



#testing joeys function with n=2000
n <- 2000
alpha <- 5
beta <- 7
lambda <- 1.5
params = c(lambda, alpha, beta) #needed for some functions, not simulation function
# 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")
sim1


#  still n = 2000
# alpha, beta, lambda already defined above
T=sim1[500]
tao = sim1[1:500]
omega=sim1[2000]
PCSend<-PCS(lambda,beta,alpha,T,tao,omega)
PCSend
#expecting m = 2000 as a predicted cascade size
# got 2208.175
# ahh not too close

#plot this
plot(sim1, y,  type="l", ylim=c(0,2200),  main="Hawkes Process Simulation with Predicted Cascade Size (n = 2000)", 
     ylab="Arrivals", xlab="Time")
PCSline = c(500,PCSend)
PCSy=c(T,omega)
lines(PCSy, PCSline, col="blue")
legend("topleft",legend=c("Original Simulation","Predicted Cascade Size (Estimated after arrival 500)"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))



#  still n = 2000
# alpha, beta, lambda already defined above
T=sim1[1000]
tao = sim1[1:1000]
omega=sim1[2000]
PCSend<-PCS(lambda,beta,alpha,T,tao,omega)
PCSend
#expecting m = 2000 as a predicted cascade size
# got 2192.629
# ahh not too close

#plot this
plot(sim1, y,  type="l", ylim=c(0,2200),  main="Hawkes Process Simulation with Predicted Cascade Size (n = 2000)", 
     ylab="Arrivals", xlab="Time")
PCSline = c(1000,PCSend)
PCSy=c(T,omega)
lines(PCSy, PCSline, col="blue")
legend("topleft",legend=c("Original Simulation","Predicted Cascade Size (Estimated after arrival 1000)"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))



#  still n = 2000
# show the function 1500 and predict the last 500
T=sim1[1500]
tao = sim1[1:1500]
omega=sim1[2000]
PCSend<-PCS(lambda,beta,alpha,T,tao,omega)
PCSend
#expecting m = 2000 as a predicted cascade size
# got 2073.75
# pretty close

#plot this
plot(sim1, y,  type="l", ylim=c(0,2200), main="Hawkes Process Simulation with Predicted Cascade Size (n = 2000)", 
     ylab="Arrivals", xlab="Time")
PCSline = c(1500,PCSend)
PCSy=c(T,omega)
lines(PCSy, PCSline, col="blue")
legend("topleft",legend=c("Original Simulation","Predicted Cascade Size (Estimated after arrival 1500)"),
       text.col=c("black","blue"),pch=c(16,15),col=c("black","blue"))