diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index 15d148d1..5de74e98 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -230,41 +230,40 @@ First, construct our $F$ from {eq}`dfcvsde` and $G$ from {eq}`dG` ```{code-cell} julia function F(x, p, t) - s, i, r, d, R_0, delta = x - (;gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p - - return [-gamma * R_0 * s * i; # ds/dt - gamma * R_0 * s * i - gamma * i; # di/dt - (1 - delta) * gamma * i; # dr/dt - delta * gamma * i; # dd/dt - eta * (R_bar_0(t, p) - R_0); # dR_0/dt - theta * (delta_bar - delta); # ddelta/dt - ] + s, i, r, d, R₀, δ = x + (;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p + + return [-γ*R₀*s*i; # ds/dt + γ*R₀*s*i - γ*i; # di/dt + (1-δ)*γ*i; # dr/dt + δ*γ*i; # dd/dt + η*(R̄₀(t, p) - R₀);# dR₀/dt + θ*(δ_bar - δ); # dδ/dt + ] end function G(x, p, t) - s, i, r, d, R_0, delta = x - (;gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p + s, i, r, d, R₀, δ = x + (;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p - return [0; 0; 0; 0; sigma * sqrt(R_0); xi * sqrt(delta * (1 - delta))] + return [0; 0; 0; 0; σ*sqrt(R₀); ξ*sqrt(δ * (1-δ))] end ``` Next create a settings generator, and then define a [SDEProblem](https://docs.sciml.ai/stable/tutorials/sde_example/#Example-2:-Systems-of-SDEs-with-Diagonal-Noise-1) with Diagonal Noise. ```{code-cell} julia -p_gen(;T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20, - R_0_n = 1.6, R_bar_0 = (t, p) -> p.R_0_n, delta_bar = 0.01, - sigma = 0.03, xi = 0.004, theta = 0.2, N = 3.3E8) = (;T, gamma, eta, R_0_n, R_bar_0, delta_bar, sigma, xi, theta, N) - -p = p_gen() # use all defaults +p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, η = 1.0 / 20, + R₀_n = 1.6, R̄₀ = (t, p) -> p.R₀_n, δ_bar = 0.01, + σ = 0.03, ξ = 0.004, θ = 0.2, N = 3.3E8) +p = p_gen() # use all defaults i_0 = 25000 / p.N r_0 = 0.0 d_0 = 0.0 s_0 = 1.0 - i_0 - r_0 - d_0 -R_bar_0_0 = 0.5 # starting in lockdown -delta_0 = p.delta_bar -x_0 = [s_0, i_0, r_0, d_0, R_bar_0_0, delta_0] +R̄₀_0 = 0.5 # starting in lockdown +δ_0 = p.δ_bar +x_0 = [s_0, i_0, r_0, d_0, R̄₀_0, δ_0] prob = SDEProblem(F, G, x_0, (0, p.T), p) ``` @@ -371,29 +370,29 @@ We will shut down the shocks to the mortality rate (i.e. $\xi = 0$) to focus on Consider $\eta = 1/50$ and $\eta = 1/20$, where we start at the same initial condition of $R_0(0) = 0.5$. ```{code-cell} julia -function generate_eta_experiment(eta; p_gen = p_gen, trajectories = 100, saveat = 1.0, x_0 = x_0, T = 120.0) - p = p_gen(eta = eta, xi = 0.0) +function generate_η_experiment(η; p_gen = p_gen, trajectories = 100, + saveat = 1.0, x_0 = x_0, T = 120.0) + p = p_gen(η = η, ξ = 0.0) ensembleprob = EnsembleProblem(SDEProblem(F, G, x_0, (0, T), p)) - sol = solve(ensembleprob, SOSRI(), EnsembleThreads(), trajectories = trajectories, saveat = saveat) + sol = solve(ensembleprob, SOSRI(), EnsembleThreads(), + trajectories = trajectories, saveat = saveat) return EnsembleSummary(sol) end # Evaluate two different lockdown scenarios -eta_1 = 1 / 50 -eta_2 = 1 / 20 -summ_1 = generate_eta_experiment(eta_1) -summ_2 = generate_eta_experiment(eta_2) - -plot(summ_1, idxs = (4, 5), +η_1 = 1/50 +η_2 = 1/20 +summ_1 = generate_η_experiment(η_1) +summ_2 = generate_η_experiment(η_2) +plot(summ_1, idxs = (4,5), title = ["Proportion Dead" L"R_0"], ylabel = [L"d(t)" L"R_0(t)"], xlabel = L"t", - legend = :topleft, - labels = L"Middle 95% Quantile, $\eta = %$eta_1$", + legend = [:topleft :bottomright], + labels = L"Middle 95% Quantile, $\eta = %$η_1$", layout = (2, 1), size = (900, 900), fillalpha = 0.5) - -plot!(summ_2, idxs = (4, 5), - legend = :topleft, - labels = L"Middle 95% Quantile, $\eta = %$eta_2$", size = (900, 900), fillalpha = 0.5) +plot!(summ_2, idxs = (4,5), + legend = [:topleft :bottomright], + labels = L"Middle 95% Quantile, $\eta = %$η_2$", size = (900, 900), fillalpha = 0.5) ``` While the the mean of the $d(t)$ increases, unsurprisingly, we see that the 95% quantile for later time periods is also much larger - even after the $R_0$ has converged. @@ -412,24 +411,24 @@ Since empirical estimates of $R_0(t)$ discussed in {cite}`NBERw27128` and other We start the model with 100,000 active infections. ```{code-cell} julia -R_0_L = 0.5 # lockdown -eta_experiment = 1.0 / 10 -sigma_experiment = 0.04 +R₀_L = 0.5 # lockdown +η_experiment = 1.0/10 +σ_experiment = 0.04 +R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0 +R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0 -R_bar_0_lift_early(t, p) = t < 30.0 ? R_0_L : 2.0 -R_bar_0_lift_late(t, p) = t < 120.0 ? R_0_L : 2.0 +p_early = p_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment) +p_late = p_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment) -p_early = p_gen(R_bar_0 = R_bar_0_lift_early, eta = eta_experiment, sigma = sigma_experiment) -p_late = p_gen(R_bar_0 = R_bar_0_lift_late, eta = eta_experiment, sigma = sigma_experiment) # initial conditions i_0 = 100000 / p_early.N r_0 = 0.0 d_0 = 0.0 s_0 = 1.0 - i_0 - r_0 - d_0 -delta_0 = p_early.delta_bar +δ_0 = p_early.δ_bar -x_0 = [s_0, i_0, r_0, d_0, R_0_L, delta_0] # start in lockdown +x_0 = [s_0, i_0, r_0, d_0, R₀_L, δ_0] # start in lockdown prob_early = SDEProblem(F, G, x_0, (0, p_early.T), p_early) prob_late = SDEProblem(F, G, x_0, (0, p_late.T), p_late) ``` @@ -439,13 +438,13 @@ Simulating for a single realization of the shocks, we see the results are qualit ```{code-cell} julia sol_early = solve(prob_early, SOSRI()) sol_late = solve(prob_late, SOSRI()) -plot(sol_early, vars = [5, 1, 2, 4], +plot(sol_early, vars = [5, 1,2,4], title = [L"R_0" "Susceptible" "Infected" "Dead"], layout = (2, 2), size = (900, 600), ylabel = [L"R_0(t)" L"s(t)" L"i(t)" L"d(t)"], xlabel = L"t", legend = [:bottomright :topright :topright :topleft], label = ["Early" "Early" "Early" "Early"]) -plot!(sol_late, vars = [5, 1, 2, 4], +plot!(sol_late, vars = [5, 1,2,4], legend = [:bottomright :topright :topright :topleft], label = ["Late" "Late" "Late" "Late"], xlabel = L"t") ``` @@ -536,26 +535,24 @@ We will redo the "Ending Lockdown" simulation from above, where the only differe ```{code-cell} julia function F_reinfect(x, p, t) - s, i, r, d, R_0, delta = x - (;gamma, R_bar_0, eta, sigma, xi, theta, delta_bar, nu) = p - - return [ - -gamma * R_0 * s * i + nu * r; # ds/dt - gamma * R_0 * s * i - gamma * i; # di/dt - (1 - delta) * gamma * i - nu * r; # dr/dt - delta * gamma * i; # dd/dt - eta * (R_bar_0(t, p) - R_0); # dR_0/dt - theta * (delta_bar - delta); # ddelta/dt - ] + s, i, r, d, R₀, δ = x + (;γ, R̄₀, η, σ, ξ, θ, δ_bar, ν) = p + + return [-γ*R₀*s*i + ν*r; # ds/dt + γ*R₀*s*i - γ*i; # di/dt + (1-δ)*γ*i - ν*r # dr/dt + δ*γ*i; # dd/dt + η*(R̄₀(t, p) - R₀);# dR₀/dt + θ*(δ_bar - δ); # dδ/dt + ] end -p_re_gen(;T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20, - R_0_n = 1.6, R_bar_0 = (t, p) -> p.R_0_n, - delta_bar = 0.01, sigma = 0.03, xi = 0.004, theta = 0.2, - N = 3.3E8, nu = 1/360) = (;T, gamma, eta, R_0_n, R_bar_0, delta_bar, sigma, xi, theta, N, nu) +p_re_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, η = 1.0 / 20, + R₀_n = 1.6, R̄₀ = (t, p) -> p.R₀_n, + δ_bar = 0.01, σ = 0.03, ξ = 0.004, θ = 0.2, N = 3.3E8, ν = 1/360) -p_re_early = p_re_gen(R_bar_0 = R_bar_0_lift_early, eta = eta_experiment, sigma = sigma_experiment) -p_re_late = p_re_gen(R_bar_0 = R_bar_0_lift_late, eta = eta_experiment, sigma = sigma_experiment) +p_re_early = p_re_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment) +p_re_late = p_re_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment) trajectories = 400 saveat = 1.0 @@ -608,4 +605,3 @@ In this case, there are significant differences between the early and late death This bleak simulation has assumed that no individuals has long-term immunity and that there will be no medical advancements on that time horizon - both of which are unlikely to be true. Nevertheless, it suggests that the timing of lifting lockdown has a more profound impact after 18 months if we allow stochastic shocks imperfect immunity. - diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index 06606dd7..021af008 100644 --- a/lectures/continuous_time/seir_model.md +++ b/lectures/continuous_time/seir_model.md @@ -171,13 +171,13 @@ We begin by implementing a simple version of this model with a constant $R_0$ an First, define the system of equations ```{code-cell} julia -function F_simple(x, p, t; gamma = 1/18, R_0 = 3.0, sigma = 1/5.2) +function F_simple(x, p, t; γ = 1/18, R₀ = 3.0, σ = 1/5.2) s, e, i, r = x - return [-gamma * R_0 * s * i; # ds/dt - gamma * R_0 * s * i - sigma * e; # de/dt - sigma * e - gamma * i; # di/dt - gamma * i; # dr/dt + return [-γ*R₀*s*i; # ds/dt = -γR₀si + γ*R₀*s*i - σ*e;# de/dt = γR₀si -σe + σ*e - γ*i; # di/dt = σe -γi + γ*i; # dr/dt = γi ] end ``` @@ -195,7 +195,7 @@ tspan = (0.0, 350.0) # ≈ 350 days prob = ODEProblem(F_simple, x_0, tspan) ``` -With this, choose an ODE algorithm and solve the initial value problem. A good default algorithm for non-stiff ODEs of this sort might be `Tsit5()`, which is the Tsitouras 5/4 Runge-Kutta method. +With this, choose an ODE algorithm and solve the initial value problem. A good default algorithm for non-stiff ODEs of this sort might be `Tsit5()`, which is the Tsitouras 5/4 Runge-Kutta method). ```{code-cell} julia sol = solve(prob, Tsit5()) @@ -242,7 +242,7 @@ While we could integrate the deaths given the solution to the model ex-post, it This is a common trick when solving systems of ODEs. While equivalent in principle to using the appropriate quadrature scheme, this becomes especially convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is not a regular time grid). Note that when doing so, $d(0) = \int_0^0 \delta \gamma i(\tau) d \tau = 0$ is the initial condition. -The system {eq}`seir_system` and the supplemental equations can be written in vector form $x := [s, e, i, r, R_0, c, d]$ with parameter tuple $p := (\sigma, \gamma, \eta, \delta, \bar{R}_0(\cdot))$ +The system {eq}`seir_system` and the supplemental equations can be written in vector form $x := [s, e, i, r, R₀, c, d]$ with parameter tuple $p := (\sigma, \gamma, \eta, \delta, \bar{R}_0(\cdot))$ Note that in those parameters, the targeted reproduction number, $\bar{R}_0(t)$, is an exogenous function. @@ -289,36 +289,35 @@ First, construct our $F$ from {eq}`dfcv` ```{code-cell} julia function F(x, p, t) - s, e, i, r, R_0, c, d = x - (;sigma, gamma, R_bar_0, eta, delta) = p - - return [-gamma * R_0 * s * i; # ds/dt - gamma * R_0 * s * i - sigma * e; # de/dt - sigma * e - gamma * i; # di/dt - gamma * i; # dr/dt - eta * (R_bar_0(t, p) - R_0); # dR_0/dt - sigma * e; # dc/dt - delta * gamma * i; # dd/dt + s, e, i, r, R₀, c, d = x + (;σ, γ, R̄₀, η, δ) = p + + return [-γ*R₀*s*i; # ds/dt + γ*R₀*s*i - σ*e; # de/dt + σ*e - γ*i; # di/dt + γ*i; # dr/dt + η*(R̄₀(t, p) - R₀);# dR₀/dt + σ*e; # dc/dt + δ*γ*i; # dd/dt ] end; - ``` This function takes the vector `x` of states in the system and extracts the fixed parameters passed into the `p` object. -The only confusing part of the notation is the `R_bar_0(t, p)` which evaluates the `p.R_bar_0` at this time (and also allows it to depend on the `p` parameter). +The only confusing part of the notation is the `R̄₀(t, p)` which evaluates the `p.R̄₀` at this time (and also allows it to depend on the `p` parameter). ### Parameters -#####HERE ENDED + The baseline parameters are put into a named tuple generator (see previous lectures using [Parameters.jl](https://github.com/mauro3/Parameters.jl)) with default values discussed above. ```{code-cell} julia -p_gen(;T = 550.0, gamma = 1.0 / 18, sigma = 1 / 5.2, eta = 1.0 / 20, - R_0_n = 1.6, delta = 0.01, N = 3.3E8, - R_bar_0 = (t, p) -> p.R_0_n) = (;T, gamma, sigma, eta, R_0_n, delta, N, R_bar_0) +p_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20, + R₀_n = 1.6, δ = 0.01, N = 3.3E8, + R̄₀ = (t, p) -> p.R₀_n); ``` -Note that the default $\bar{R}_0(t)$ function always equals $R_{0n}$ -- a parameterizable natural level of $R_0$ used only by the `R_bar_0` function +Note that the default $\bar{R}_0(t)$ function always equals $R_{0n}$ -- a parameterizable natural level of $R_0$ used only by the `R̄₀` function Setting initial conditions, we choose a fixed $s, i, e, r$, as well as $R_0(0) = R_{0n}$ and $m(0) = 0.01$ @@ -329,7 +328,7 @@ i_0 = 1E-7 e_0 = 4.0 * i_0 s_0 = 1.0 - i_0 - e_0 -x_0 = [s_0, e_0, i_0, 0.0, p.R_0_n, 0.0, 0.0] +x_0 = [s_0, e_0, i_0, 0.0, p.R₀_n, 0.0, 0.0] tspan = (0.0, p.T) prob = ODEProblem(F, x_0, tspan, p) ``` @@ -373,9 +372,9 @@ Let's start with the case where $\bar{R}_0(t) = R_{0n}$ is constant. We calculate the time path of infected people under different assumptions of $R_{0n}$: ```{code-cell} julia -R_0_n_vals = range(1.6, 3.0, length = 6) -sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R_0_n = R_0_n)), - Tsit5(), saveat=0.5) for R_0_n in R_0_n_vals]; +R₀_n_vals = range(1.6, 3.0, length = 6) +sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R₀_n = R₀_n)), + Tsit5(), saveat=0.5) for R₀_n in R₀_n_vals]; ``` Here we chose `saveat=0.5` to get solutions that were evenly spaced every `0.5`. @@ -385,7 +384,7 @@ Changing the saved points is just a question of storage/interpolation, and does Let's plot current cases as a fraction of the population. ```{code-cell} julia -labels = permutedims([L"R_0 = %$r" for r in R_0_n_vals]) +labels = permutedims([L"R_0 = %$r" for r in R₀_n_vals]) infecteds = [sol[3,:] for sol in sols] plot(infecteds, label=labels, legend=:topleft, lw = 2, xlabel = L"t", ylabel = L"i(t)", title = "Current Cases") @@ -414,21 +413,21 @@ In the simple case, where $\bar{R}_0(t) = R_{0n}$ is independent of the state, t We will examine the case where $R_0(0) = 3$ and then it falls to $R_{0n} = 1.6$ due to the progressive adoption of stricter mitigation measures. -The parameter `eta` controls the rate, or the speed at which restrictions are +The parameter `η` controls the rate, or the speed at which restrictions are imposed. We consider several different rates: ```{code-cell} julia -eta_vals = [1/5, 1/10, 1/20, 1/50, 1/100] -labels = permutedims([L"\eta = %$eta" for eta in eta_vals]); +η_vals = [1/5, 1/10, 1/20, 1/50, 1/100] +labels = permutedims([L"\eta = %$η" for η in η_vals]); ``` Let's calculate the time path of infected people, current cases, and mortality ```{code-cell} julia x_0 = [s_0, e_0, i_0, 0.0, 3.0, 0.0, 0.0] -sols = [solve(ODEProblem(F, x_0, tspan, p_gen(eta=eta)), Tsit5(), saveat=0.5) for eta in eta_vals]; +sols = [solve(ODEProblem(F, x_0, tspan, p_gen(η=η)), Tsit5(), saveat=0.5) for η in η_vals]; ``` Next, plot the $R_0$ over time: @@ -469,11 +468,11 @@ The parameters considered here start the model with 25,000 active infections and 75,000 agents already exposed to the virus and thus soon to be contagious. ```{code-cell} julia -R_0_L = 0.5 # lockdown -R_bar_0_lift_early(t, p) = t < 30.0 ? R_0_L : 2.0 -R_bar_0_lift_late(t, p) = t < 120.0 ? R_0_L : 2.0 -p_early = p_gen(R_bar_0 = R_bar_0_lift_early, eta = 10.0) -p_late = p_gen(R_bar_0 = R_bar_0_lift_late, eta = 10.0) +R₀_L = 0.5 # lockdown +R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0 +R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0 +p_early = p_gen(R̄₀= R̄₀_lift_early, η = 10.0) +p_late = p_gen(R̄₀= R̄₀_lift_late, η = 10.0) # initial conditions @@ -481,9 +480,9 @@ i_0 = 25000 / p_early.N e_0 = 75000 / p_early.N s_0 = 1.0 - i_0 - e_0 -x_0 = [s_0, e_0, i_0, 0.0, R_0_L, 0.0, 0.0] # start in lockdown +x_0 = [s_0, e_0, i_0, 0.0, R₀_L, 0.0, 0.0] # start in lockdown -# create two problems, with rapid movement of R_0(t) towards R_bar_0(t) +# create two problems, with rapid movement of R₀(t) towards R̄₀(t) prob_early = ODEProblem(F, x_0, tspan, p_early) prob_late = ODEProblem(F, x_0, tspan, p_late) ``` @@ -502,7 +501,7 @@ plot!(sol_late, vars = [7], label = "Lift Late", xlabel = L"t") Next we examine the daily deaths, $\frac{d D(t)}{dt} = N \delta \gamma i(t)$. ```{code-cell} julia -flow_deaths(sol, p) = p.N * p.delta * p.gamma * sol[3,:] +flow_deaths(sol, p) = p.N * p.δ * p.γ * sol[3,:] plot(sol_early.t, flow_deaths(sol_early, p_early), title = "Flow Deaths", label = "Lift Early") plot!(sol_late.t, flow_deaths(sol_late, p_late), label = "Lift Late", xlabel = L"t") @@ -518,4 +517,3 @@ Despite its richness, the model above is fully deterministic. The policy $\bar{ One way that randomness can lead to aggregate fluctuations is the granularity that comes through the discreteness of individuals. This topic, the connection between SDEs and the Langevin equations typically used in the approximation of chemical reactions in well-mixed media is explored in further lectures on continuous time Markov chains. Instead, in the {doc}`next lecture `, we will concentrate on randomness that comes from aggregate changes in behavior or policy. -