### read data

global_case <- read.csv("/YOUR PATH/COVID_global_case.csv",header=T)
global_death <- read.csv("/YOUR PATH/COVID_global_death.csv",header=T)
global_recover <- read.csv("/YOUR PATH/COVID_global_recover.csv",header=T)

# Hubei

temp <- ncol(subset(global_case, Province.State=="Hubei"))

case_hubei <- as.numeric(subset(global_case, Province.State=="Hubei")[5:temp])
death_hubei <- as.numeric(subset(global_death, Province.State=="Hubei")[5:temp])
recover_hubei <- as.numeric(subset(global_recover, Province.State=="Hubei")[5:temp])

CFR_hubei <- death_hubei/(death_hubei+recover_hubei)

# China outside Hubei

case_china_matrix <- subset(global_case, Country.Region=="China")
death_china_matrix <- subset(global_death, Country.Region=="China")
recover_china_matrix <- subset(global_recover, Country.Region=="China")

case_china <- colSums(case_china_matrix[,5:(ncol(case_china_matrix))])
death_china <- colSums(death_china_matrix[,5:(ncol(case_china_matrix))])
recover_china <- colSums(recover_china_matrix[,5:(ncol(case_china_matrix))])

case_outside_hubei <- as.numeric(case_china-case_hubei)
death_outside_hubei <- as.numeric(death_china-death_hubei)
recover_outside_hubei <- as.numeric(recover_china-recover_hubei)

CFR_outside_hubei <- death_outside_hubei/(death_outside_hubei+recover_outside_hubei)

CFR_estimate_temp <- rep(0, 80) # store CFR for different Y

for (Y in 5:80){ ### CFR estimate for different Y

###### observed data

case <- case_hubei[1:Y]
death <- death_hubei[1:Y]
recover <- recover_hubei[1:Y]

### compute new case, new death, and new recover

new_case <- case
new_death <- death
new_recover <- recover

for (i in 2:length(case)){
	new_case[i] <- case[i] - case[i-1]
	new_death[i] <- death[i] - death[i-1]
	new_recover[i] <- recover[i] - recover[i-1]
}

### start function

function_Poisson_CFR <- function(CFR, lambda_death, lambda_recover){

### create death matrix

death_matrix <- matrix(0, nrow=length(case)-1, ncol=length(case)-1)

for (i in 1:nrow(death_matrix)){
	for (j in 1:i){
		death_matrix[i,j] <- CFR*dpois(j-1,lambda_death)*new_case[i-j+1]
	}
}

### create recover matrix

recover_matrix <- matrix(0, nrow=length(case)-1, ncol=length(case)-1)

for (i in 1:nrow(recover_matrix)){
	for (j in 1:i){
		recover_matrix[i,j] <- (1-CFR)*dpois(j-1,lambda_recover)*new_case[i-j+1]
	}
}

### compute expected death and recover

expect_death <- rowSums(death_matrix)
expect_recover <- rowSums(recover_matrix)

### compute overall chi-square

chisq <- 0

for (i in 2:length(case)){
	chisq <- chisq + (new_death[i]-expect_death[i-1])^2/expect_death[i-1]
	chisq <- chisq + (new_recover[i]-expect_recover[i-1])^2/expect_recover[i-1]
}

return(chisq) ### return the overall chi-square

} ### end function

### write function_Poisson_CFR as a 1-argument function

f_test <- function(x) function_Poisson_CFR(x[1],x[2],x[3])

### optimization

test_output <- optim(par=c(0.01,4,11),fn=f_test,lower=c(0.0001,1,5),upper=c(0.9999,15,30),method = "L-BFGS-B")

CFR_estimate_temp[Y] <- test_output$par[1]

} ### CFR for different Y ends

CFR_estimate_temp