Skip to content

Commit

Permalink
update so will merge
Browse files Browse the repository at this point in the history 
***ORDER matters make sure to merge continuous time last
  • Loading branch information
@annabellasd
annabellasd committed Oct 11, 2023
1 parent befd388 commit c1f5ed2
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 107 deletions.
126 changes: 61 additions & 65 deletions lectures/continuous_time/covid_sde.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```
Expand Down Expand Up @@ -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.
Expand All @@ -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)
```
Expand All @@ -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")
```
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.

82 changes: 40 additions & 42 deletions lectures/continuous_time/seir_model.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
Expand All @@ -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())
Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -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$

Expand All @@ -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)
```
Expand Down Expand Up @@ -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`.
Expand All @@ -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")
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -469,21 +468,21 @@ 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
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)
```
Expand All @@ -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")
Expand All @@ -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 <covid_sde>`, we will concentrate on randomness that comes from aggregate changes in behavior or policy.

0 comments on commit c1f5ed2

Please sign in to comment.