# go back to Flu Mortality

#Excess deaths

diff = c(1.4486409, 9.2955840, 0.8760308, 14.6437410, 28.4196368, 47.2210290, 73.0648590, 136.9674535, 166.9442884, 163.0097641, 125.1769971, 70.4576305, 44.8616645, 71.3973124, 34.0708812, 76.8866802, 36.8469589, 36.9518742, 22.1994882, 27.5857962, 4.1047848, 11.7485187, 10.5072564, 12.3695902, 25.3226122, 15.3521012, 8.5782790, 1.7418885, 10.9162854, 16.3306561, 24.3761382, 3.2338755, 10.0183877, 9.6904685, 2.2399656, 6.6585022, 2.9395978, 9.6335156)

#Time in weeks of excess deaths

tb = c(0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 28, 29, 30, 33, 34, 36, 37, 38, 39, 40, 41, 47)

#SIR model

SIR = function(Time, State, Pars) {

with(as.list(c(State, Pars)), {
	
	DecreaseS = -1 * beta * S * I
	DecreaseI = -1 * alpha * I
	
	dS = DecreaseS
	dI = -1 * DecreaseS + DecreaseI
	
	return(list(c(dS, dI)))
	})

}

#load the deSolve package

#Log-likelihood function

logl = function(beta1, beta2, I10, I20) {

time = seq(0, 329, 1)

pars1 = c(beta = beta1, alpha = (1/7))
state1 = c(S = (27629291 - I10), I = I10)
out1 = ode(y = state1, times = time, func = SIR, parms = 	pars1)
S1 = out1[,"S"]

pars2 = c(beta = beta2, alpha = (1/7))
state2 = c(S = (27629291 - I20), I = I20)
out2 = ode(y = state2, times = time, func = SIR, parms = 	pars2)
S2 = out2[,"S"]

lamb1f = function(t)
{
	if(t == 0)
	{
		S1val = S1[t]
		
		return(S1val)
	}
	else
	{
		S1val = S1[t*7] - S1[7*(t - 1) + 1]
	
		return(S1val)
	}
}

lamb2f = function(t)
{
	if(t == 0)
	{
		S2val = S2[t]
		
		return(S2val)
	}
	else
	{
		S2val = S2[t*7] - S2[7*(t - 1) + 1]
	
		return(S2val)
	}
}

lamb1 = c()
lamb2 = c()

for(i in 1:38)
{
	lamb1[i] = lamb1f(tb[i])
}

for(i in 1:38)
{
	lamb2[i] = lamb2f(tb[i])
}

loglik = sum(diff*log(lamb1 + lamb2, base = exp(1)) - lamb1 - lamb2)

return(loglik)

}

#Optimization solver

optim(c(0.258, 0.236, 80000, 30000), fn = logl, method=“L-BFGS-B”, lower = c(0, 0, 1, 1), upper = c(Inf, Inf, 27629291, 27629291))