Quick implementation of SEIHR in Python (no differentiation between households and other yet):
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
N = 6.6e6 #population size
N_INFECTED = N/1e3
def covid19(x,t):
a1 = 0.1/14*100 #probability of contracting from infected (10% risk of infection if in contact for 14 days, 100 contacts within this period)
a2 = a1/10 #probability of contracting from presymptomatic (during virus incubation period)
d = 1/21 #rate of recovery from infection without hospitalisation (average recovery in 2.5 weeks)
e = 1/35 #rate of complete recovery at hospital
f = e*0.16 #rate of death at hospital (only hospitalised die in this model)
g = d*0.2 #rate of intected degrading to critical cases requiring hospitalisation (20% requires hospitalisation)
k = 1/4 #rate of developing infection, where the average incubation time is 1/k
S = x[0]
E = x[1]
I = x[2]
H = x[3]
R = x[4]
dSdt = -a1/N*S*I - a2/N*S*E #susceptible stock
dEdt = a1/N*S*I + a2/N*S*E - k*E #exposed stock
dIdt = k*E - g*I - d*I #infected stock
dHdt = g*I - e*H - f*H #hospitalised stock
dRdt = d*I + e*H #recovered stock
return [dSdt,dEdt,dIdt,dHdt,dRdt]
x0 = [N-3*N_INFECTED,N_INFECTED,N_INFECTED,N_INFECTED,0]
t0 = 0
t1 = 180
dt = 0.01
t = np.linspace(t0,t1,int((t1-t0)/dt))
x = odeint(covid19,x0,t)
S = x[:,0]
E = x[:,1]
I = x[:,2]
H = x[:,3]
R = x[:,4]
D = N-(S+E+I+R+H)
plt.figure(figsize = (10,10))
plt.plot(t,S,label='Susceptible')
plt.plot(t,E,label='Exposed')
plt.plot(t,I,label='Infected')
plt.plot(t,R,label='Recovered')
plt.plot(t,H,label='Hospitalised')
plt.plot(t,D, label='Deceased')
plt.legend()
plt.figure(figsize = (10,10))
plt.plot(t,D/(R+D),label='Fatality rate')