## com55a: a more polished version of com117* from nilswork22

#  i run this code for two occasions:
#  Nor-Den 27-25 nov 2022
#  Nor-Den 23-23 dec 2023 

# com117c of december 2023,
# same method as Real Time Real Excitement Plots on earlier occasions 

# Denmark Norway women's handball,
# how exciting was it, really? 23 23 after 60 minutes 

# here plain independent Poissons  

# can I also do the more complex model with (lambda, phi)?

# time is in minutes 
# goals per team in 60 minutes: about 27 

lambda <- 27/60 
# gamma <- 0.418 
# tmax <- 60 

## what's the probability of a draw? 

xval <- seq(0,100) 
checkval <- dpois(xval,60*lambda)^2 

plot(xval,checkval) 
probdraw <- sum(checkval) # is prob(X(60) = Y(60)) is 0.0544 

# so 60*lambda*2 = 54 expected goals in a match  
# and 60*lambda = 27 is Poisson parameter per team 

# handball <- matrix(scan("denmark_norway_handballdec2023",skip=6),byrow=T,ncol=4)
handball <- matrix(scan("denmark_norway_handballnov2022",skip=6),byrow=T,ncol=4)

time <- handball[ ,1] + handball[ ,2]/60 # times for goals 
den <- handball[ ,3]
nor <- handball[ ,4]

zz <- den + nor 
cbind(handball,den,nor,zz)

## ok, Norway is A, Denmark is B 

nnA <- max(nor) # 23 
nnB <- max(den) # 23 

nn <- nnA + nnB # 46 goals at time points given in time 

timeplus <- c(0,time,60) # full window, adding on 0.00 and 60.00
nnplus <- nn + 2 

teamA0 <- nor # length nn
teamB0 <- den # length nn  

teamA <- c(0,teamA0,nnA)
teamB <- c(0,teamB0,nnB)

cbind(timeplus,teamA,teamB) # starting at 0-0 after t = 0 

# testing that Poisson rate is the same: 

## making handball_finale_fig12.pdf : 

par(bty="l")
par(cex.lab=1.4142) 

matplot(timeplus,cbind(teamA,teamB),type="l",lwd=3,
   xlab="match time", ylab="Norway and Denmark")

matplot(timeplus,teamB/25,type="l")

## keep <- cbind(time,den,nor)

par(bty="l")
par(cex.lab=1.4142) 

matplot(timeplus,cbind(teamA,teamB),type="l",
   xlab="minutes",ylab="Norway vs Denmark", lwd=3)

matplot(timeplus,teamA - teamB,type="l")
matlines(timeplus,0*timeplus,type="l",lty=2)

probA <- 0*timeplus # nor
probB <- 0*timeplus # den 
probD <- 0*timeplus # draw 

# at time zero: 
probD[1] <- probdraw
probA[1] <- 0.5*(1-probdraw)
probB[1] <- 0.5*(1-probdraw)

# for final figure: sim = 10^6 
sim <- 10^6

for (k in 2:(nn+1))
{
tnow <- timeplus[k]
thresholdnow <- teamB[k] - teamA[k] 
# 
lambdanow <- (60-tnow)*lambda 
# 
# simulating A' - B' 
diffsim <- rpois(sim,lambdanow) - rpois(sim,lambdanow)
# 
probA[k] <- mean(1*(diffsim > thresholdnow))
probB[k] <- mean(1*(diffsim < thresholdnow))
probD[k] <- mean(1*(diffsim == thresholdnow))
# 
print(c(k,probA[k]))
}

probA[nn+2] <- probA[nn+1]
probB[nn+2] <- probB[nn+1]
probD[nn+2] <- probD[nn+1]

# some polish for the case of nov 2022, where Nor-Den 27-25: 
probA[nn+2] <- 1
probB[nn+2] <- 0
probD[nn+2] <- 0 

# some polish for the case of dec 2023, where there is a draw: 
# probA[nn+2] <- 0
# probB[nn+2] <- 0
# probD[nn+2] <- 1 

## making handball_finale_fig11.pdf : 

par(bty="l")
par(cex.lab=1.4142)

matplot(timeplus,cbind(probA,probB,probD),type="l",
   lwd=3,lty=c(1,2,3),xlab="minutes",
   ylab="Pr(Nor wins), Pr(Den wins), Pr(draw)")

matlines(timeplus,0*timeplus+0.5,type="l",lty=3,lwd=3)

## making handball_finale_fig11B22.pdf : one for November 2022 
## making handball_finale_fig11B23.pdf : one for December 2023 

par(bty="l")
par(cex.lab=1.4142)

matplot(timeplus,cbind(probA,probB,probD),type="l",
   lwd=3,lty=c(1,1,2),col=c(1,2,4),xlab="minutes",
   ylab="Pr(Nor wins), Pr(Den wins), Pr(draw)")
grid(lwd=3) 

matlines(timeplus,0*timeplus+0.5,type="l",lty=3,lwd=3)

## A alone and B alone: 

checkA <- 0*(1:nn) # 52 goals: at which of these times did A score? 
checkB <- 0*(1:nn) # 52 goals: at which of these times did B score? 

for (j in 2:nn) 
{
checkA[j] <- 1*(teamA0[j] > teamA0[j-1])
checkB[j] <- 1*(teamB0[j] > teamB0[j-1])
}

checkA[1] <- 1*(teamA0[1] == 1)
checkB[1] <- 1*(teamB0[1] == 1)

goalsA <- time[checkA == 1] # nnA times
goalsB <- time[checkB == 1] # nnB times 

tval <- seq(0,60,by=0.01)
FAval <- 0*tval
FBval <- 0*tval 

for (j in 1:length(tval))
{
FAval[j] <- mean(1*(goalsA <= tval[j]))
FBval[j] <- mean(1*(goalsB <= tval[j]))
}

matplot(tval,cbind(FAval,FBval),type="l")

ZAval <- sqrt(nnA)*(FAval - tval/60)
ZBval <- sqrt(nnB)*(FBval - tval/60)

## making handball_finale_fig13.pdf : 

par(bty="l")
par(cex.lab=1.4142) 

matplot(tval,cbind(ZAval,ZBval),type="l",lwd=3,
   xlab="match time",ylab="bridge plots, Norway and Denmark")
grid(lwd=3)