CoCalc Public FilesPublic / Covid19EAF200.ipynbOpen with one click!
Author: Caio Machado
Views : 23

Material preparado por Caio Machado para el curso Aplicaciones Matemáticas para Economía y Negocios (solamente para fines didácticos)

In [1]:
using DifferentialEquations, Parameters, Plots; # Packages que vamos usar #plotly(); # Uncomment para graficos interactivos

1. Modelo base

El modelo presentado es una versión del modelo en Fernandez-Villaverde and Jones (2020) (FJ en lo que sigue). La población total es de NN personas, y cada persona puede estar en uno de cinco estados diferentes: susceptible, infectada, en recuperación, fallecida y recuperada (las personas se mueven de un estado a otro en este orden, excepto por los últimos dos estados). Tenemos que:

St+It+Rt+Dt+Ct=N S_{t}+I_{t}+R_{t}+D_{t}+C_{t}=N

donde

St=total SusceptiblesS_{t}=\text{total Susceptibles}

It=total InfectadosI_{t}=\text{total Infectados}

Rt=total en RecuperacioˊnR_{t}=\text{total en Recuperación}

Ft=total FallecidosF_{t}=\text{total Fallecidos}

Ct=total reCuperadosC_{t}=\text{total reCuperados}

Cada período tt representa un día. Usando la notación que Δxt+1=xt+1xt\Delta x_{t+1} = x_{t+1}-x_t, la evolución del sistema es dada por:


ΔSt+1=βItSt/Nnuevas infecciones(1)\Delta S_{t+1}=-\underbrace{\beta I_{t} S_{t} / N}_{\text {nuevas infecciones}} \tag{1}
ΔIt+1=βItSt/Nnuevas infeccionesγItnuevos en recuperacioˊn(2)\Delta I_{t+1}=\underbrace{\beta I_{t} S_{t} / N}_{\text {nuevas infecciones}}-\underbrace{\gamma I_{t}}_{\text {nuevos en recuperación}} \tag{2}
ΔRt+1=γItnuevos en recuperacioˊnθRtnuevos fallecidos o recuperados(3)\Delta R_{t+1}=\underbrace{\gamma I_{t}}_{\text {nuevos en recuperación}}-\underbrace{\theta R_{t}}_{\text {nuevos fallecidos o recuperados}} \tag{3}
ΔFt+1=δθRtnuevos fallecidos(4)\Delta F_{t+1}=\underbrace{\delta \theta R_{t}}_{\text {nuevos fallecidos}} \tag{4}
ΔCt+1=(1δ)θRtnuevos recuperados(5)\Delta C_{t+1}=\underbrace{(1-\delta) \theta R_{t}}_{\text {nuevos recuperados}} \tag{5}

La interpretación de los parámetros y ecuaciones es la siguiente:

  • Todo día, cada persona infectada (ItI_t) interactúa (por un largo tiempo) en promedio con β\beta personas. Una fracción de St/NS_t/N de estas personas son susceptibles a infectarse, y por lo tanto el número total de nuevos infectados es dado por βItSt/N\beta I_{t} S_{t} / N. Por ejemplo, si el aislamiento social sube, esto seria representado por un menor valor de β\beta. Todavía, por ahora vamos a suponer que β\beta es constante en el tiempo.

  • En cada fecha, una proporción γ\gamma de los infectados se mueve al estado en recuperación. Podemos interpretar irse del estado infectado al estado en recuperación a entrar en aislamiento total, una vez que uno descubre que tiene Covid. Esta diferenciación entre infectados y en recuperación puede representar la idea que personas que tienen Covid y saben que tienen Covid infectan menos que personas que no saben que tienen Covid. En particular, estamos suponiendo que personas que saben que tienen el virus se ponen en aislamiento total y no infectan a nadie.

  • En cada fecha, una proporción θ\theta de las personas que están en recuperación dejan este estado, sea por que se han recuperado, sea por que han fallecido.

  • Una fracción δ\delta de las personas que tienen Covid fallecen, y una fracción (1δ)(1-\delta) se recuperan.

Percibe que podemos reescribir este sistema como hemos hecho en clase:

St+1=StβItSt/NS_{t+1}=S_t -\beta I_{t} S_{t} / N

It+1=It+βItSt/NγItI_{t+1}=I_t + {\beta I_{t} S_{t} / N}-{\gamma I_{t}}

Rt+1=Rt+γItθRtR_{t+1}=R_t + {\gamma I_{t}}-{\theta R_{t}}

Ft+1=Ft+δθRtF_{t+1}=F_t + {\delta \theta R_{t}}

Ct+1=Ct+(1δ)θRtC_{t+1}=C_t + {(1-\delta) \theta R_{t}}

Todavía, Julia nos pide que escribamos el sistema en el formato (1)-(5), pero es el mismo sistema por supuesto. Ahora vamos a representar este sistema dinámico utilizando Julia. Primero, definimos los parámetros del modelo:

In [2]:
@with_kw struct Para beta::Float64 gamma::Float64 delta::Float64 theta::Float64 N::Int64 end;

Los valores específicos de los parámetros vamos a determinar en cada simulación que hagamos más adelante. Ahora definimos las ecuaciones (1)-(5):

In [3]:
function f(du,u,p,t) S,I,R,F,C = u du[1] = -p.beta*I*S/p.N # Ecuación (1) du[2] = p.beta*I*S/p.N - p.gamma*I # Ecuación (2) du[3] = p.gamma*I - p.theta*R # Ecuación (3) du[4] = p.delta*p.theta*R # Ecuación (4) du[5] = (1-p.delta)*p.theta*R # Ecuación (5) end;

En todas simulaciones vamos a simular un período de 300300 días (10 meses) y además vamos a suponer que en la fecha inicial había 100 personas infectadas y los demás susceptibles (el caso 100 confirmado en Chile ocurrió en mitad de marzo). Representamos la situación inicial por el vector u0=(S0,I0,R0,F0,C0)u_0 = (S_0,I_0,R_0,F_0, C_0) abajo. El intervalo de tiempo es representado por tspantspan. Además, vamos siempre usar N=19N=19, para representar la población de Chile de aproximadamente 19 millones de personas.

In [4]:
u0 = [19-(100/1e6);100/1e6;0;0;0] # N está representado en millones, por esto dividimos por 1e6 tspan = (0.0,300.0);

Vamos a armar el problema utilizando la toolbox de Julia DifferentialEquations.jl (aunque estamos resolviendo una ecuación en diferencias, podemos usar esta toolbox para resolver problemas en tiempo discreto). Este se puede armar con 4 elementos: (i) Una ley de movimiento (la función ff); (ii) Un vector de condiciones iniciales u0u_0; (iii) Un intervalo de tiempo para las simulaciones (tspantspan); (iv) parámetros que afectan la función ff (parpar). Para armar el problema escribimos prob = ODEProblem(f,u0,tspan,par). Para resolver el problema escribimos sol = solve(prob,Euler(),dt=1), y así la variable sol almacena varias informaciones sobre la solución.

2. Simulación del modelo base

Primeramente, precisamos definir los parámetros. No voy a hacer nada muy sofisticado, la idea es solamente elegir algunos parámetros que mínimamente hagan algún sentido. Sigue como he elegido los parámetros:

  • γ\gamma: Vamos a seguir FJ y usar γ=0.2\gamma=0.2. Esto implica que en promedio, cada persona con Covid permanece en el estado infectada por 5 días (hasta pasar para el estado en recuperación).

  • θ\theta: Vamos a seguir FJ y usar θ=0.1\theta=0.1. Esto implica que cada persona infectada permanece en recuperación por 10 días, en promedio. Por lo tanto, una persona que tiene Covid lleva en promedio 15 días con el virus en nuestro modelo (5 infectada y 10 en recuperación).

  • δ\delta: La tasa de mortalidad en Chile está cerca de 1%. Es decir, el número total de fallecidos dividido por el número de casos confirmados está cerca de 1%. Todavía, es probable que haya mucha subnotificación de casos, y algunos expertos (en mi Twitter =D) estiman que el número total de casos es 5 veces mayor que el número de casos confirmados. Así que voy a usar δ=0.002\delta=0.002

  • β\beta: Por supuesto que β\beta varia en el tiempo. Por ejemplo, medidas de cuarentena obligatoria deberían reducir β\beta. Todavía, estamos suponiendo un β\beta constante (vamos a relajar este supuesto más adelante). Podemos tener una idea del valor de β\beta mirando datos que dicen cuantas personas los primeros infectados han infectado. Esto es llamado de R0\mathcal{R}_0 y es dado por:

    R0=βcontactos por dıˊa×(1/γ)promedio dıˊas infectado\mathcal{R}_0 = \underbrace{\beta}_{\text{contactos por día}} \times \underbrace{(1/\gamma)}_{\text{promedio días infectado}}

    Hay alguna variabilidad en diferentes estimaciones de R0\mathcal{R}_0, pero estas sugieren un número cerca de 2. Así usamos R0=2\mathcal{R}_0=2 que nos da β=γR0=0.4\beta=\gamma \mathcal{R}_0=0.4.

Ahora definimos los valores de los parámetros en Julia para el modelo base:

In [5]:
par = Para(beta = 0.4, gamma = 0.2, delta = 0.002, theta = 0.1, N = 19);

Finalmente podemos armar el problema dado estos parámetros y resolverlo:

In [6]:
prob1 = ODEProblem(f,u0,tspan,par) sol1 = solve(prob1,Euler(),dt=1);

Y también graficar la solución para los primeros 5 meses:

In [7]:
tplot = (0.0,150.0) # Intervalo de tiempo graficado p1 = plot(sol1,vars=1,tspan = tplot,title="Susceptibles (S_t)") p2 = plot(sol1,vars=2,tspan = tplot,title="Infectados (I_t)") p3 = plot(sol1,vars=3,tspan = tplot,title="En Recuperación (R_t)") p4 = plot(sol1,vars=4,tspan = tplot,title="Fallecidos (F_t)") p5 = plot(sol1,vars=5,tspan = tplot,title="Recuperados (C_t)") plot(plot(p2,p4,layout=(2,1)),plot(p1,p3,p5,layout=(3,1)),legend=false,xaxis="")

3. Aislamiento social

Suponga que β\beta reduce 25%. Esto representa cada persona interactuar con 75% de las personas que interactuaba antes. ¿Qué pasa con la cantidad de fallecidos y la duración de la pandemia? Para contestar esta pregunta, todo que tenemos que hacer es resolver el mismo modelo pero con otro valor de β\beta.

In [8]:
par2 = reconstruct(par, beta=0.3) # Cambiar parámetros prob2 = ODEProblem(f,u0,tspan,par2) sol2 = solve(prob2,Euler(),dt=1);

Ahora vamos graficar la solución con distintos valores de β\beta (la línea azul es el β\beta usado en el modelo base) para los primeros 8 meses:

In [9]:
tplot = (0.0,240.0) # Intervalo de tiempo graficado p1 = plot(sol1,vars=1,tspan=tplot,title="Susceptibles (S_t)") plot!(sol2,vars=1,tspan=tplot) p2 = plot(sol1,vars=2,tspan=tplot,title="Infectados (I_t)") plot!(sol2,vars=2,tspan=tplot) p3 = plot(sol1,vars=3,tspan=tplot,title="En Recuperación (R_t)") plot!(sol2,vars=3,tspan=tplot) p4 = plot(sol1,vars=4,tspan=tplot,title="Fallecidos (F_t)") plot!(sol2,vars=4,tspan=tplot) p5 = plot(sol1,vars=5,tspan=tplot,title="Recuperados (C_t)") plot!(sol2,vars=5,tspan=tplot) plot(plot(p2,p4,layout=(2,1)),plot(p1,p3,p5,layout=(3,1)),legend=false,xaxis="")