Modeling in Continuous Time Lectures

lectures/continuous_time/covid_sde.md

@@ -230,40 +230,41 @@ 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₀, δ = 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
-            ]
+    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
+           ]
 end
 
 function G(x, p, t)
-    s, i, r, d, R₀, δ = x
-    (;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p
+    s, i, r, d, R_0, delta = x
+    (; gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p
 
-    return [0; 0; 0; 0; σ*sqrt(R₀); ξ*sqrt(δ * (1-δ))]
+    return [0; 0; 0; 0; sigma * sqrt(R_0); xi * sqrt(delta * (1 - delta))]
 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 = @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
+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
 i_0 = 25000 / p.N
 r_0 = 0.0
 d_0 = 0.0
 s_0 = 1.0 - i_0 - r_0 - d_0
-R̄₀_0 = 0.5  # starting in lockdown
-δ_0 = p.δ_bar
-x_0 = [s_0, i_0, r_0, d_0, R̄₀_0, δ_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]
 
 prob = SDEProblem(F, G, x_0, (0, p.T), p)
 ```
@@ -370,29 +371,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_η_experiment(η; p_gen = p_gen, trajectories = 100,
-                              saveat = 1.0, x_0 = x_0, T = 120.0)
-    p = p_gen(η = η, ξ = 0.0)
+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)
     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
-η_1 = 1/50
-η_2 = 1/20
-summ_1 = generate_η_experiment(η_1)
-summ_2 = generate_η_experiment(η_2)
-plot(summ_1, idxs = (4,5),
+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),
     title = ["Proportion Dead" L"R_0"],
     ylabel = [L"d(t)" L"R_0(t)"], xlabel = L"t",
-    legend = [:topleft :bottomright],
-    labels = L"Middle 95% Quantile, $\eta = %$η_1$",
+    legend = :topleft,
+    labels = L"Middle 95% Quantile, $\eta = %$eta_1$",
     layout = (2, 1), 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)
+
+plot!(summ_2, idxs = (4, 5),
+    legend = :topleft,
+    labels = L"Middle 95% Quantile, $\eta = %$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.
@@ -411,24 +412,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₀_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_0_L = 0.5  # lockdown
+eta_experiment = 1.0 / 10
+sigma_experiment = 0.04
 
-p_early = p_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment)
-p_late = p_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment)
+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 = 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
-δ_0 = p_early.δ_bar
+delta_0 = p_early.delta_bar
 
-x_0 = [s_0, i_0, r_0, d_0, R₀_L, δ_0]  # start in lockdown
+x_0 = [s_0, i_0, r_0, d_0, R_0_L, delta_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)
 ```
@@ -438,13 +439,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")
 ```
@@ -535,24 +536,26 @@ 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₀, δ = 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
-            ]
+    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
+    ]
 end
 
-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_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_early = p_re_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment)
-p_re_late = p_re_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment)
+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)
 
 trajectories = 400
 saveat = 1.0

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; γ = 1/18, R₀ = 3.0, σ = 1/5.2)
+function F_simple(x, p, t; gamma = 1/18, R_0 = 3.0, sigma = 1/5.2)
     s, e, i, r = x
 
-    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
+    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
             ]
 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₀, 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_0, 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,35 +289,36 @@ First, construct our $F$ from {eq}`dfcv`
 
 ```{code-cell} julia
 function F(x, p, t)
-    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
+    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
             ]
 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̄₀(t, p)` which evaluates the `p.R̄₀` at this time (and also allows it to depend on the `p` parameter).
+The only confusing part of the notation is the `R_0(t, p)` which evaluates the `p.R_bar_0` 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 = @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);
+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)
 ```
 
-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
+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_hat_0` 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$
 
@@ -328,7 +329,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₀_n, 0.0, 0.0]
+x_0 = [s_0, e_0, i_0, 0.0, p.R_0_n, 0.0, 0.0]
 tspan = (0.0, p.T)
 prob = ODEProblem(F, x_0, tspan, p)
 ```
@@ -372,9 +373,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₀_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];
+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];
 ```
 
 Here we chose `saveat=0.5` to get solutions that were evenly spaced every `0.5`.
@@ -384,7 +385,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₀_n_vals])
+labels = permutedims([L"R_0 = %$r" for r in R_0_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")
@@ -413,21 +414,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 `η` controls the rate, or the speed at which restrictions are
+The parameter `eta` controls the rate, or the speed at which restrictions are
 imposed.
 
 We consider several different rates:
 
 ```{code-cell} julia
-η_vals = [1/5, 1/10, 1/20, 1/50, 1/100]
-labels = permutedims([L"\eta = %$η" for η in η_vals]);
+eta_vals = [1/5, 1/10, 1/20, 1/50, 1/100]
+labels = permutedims([L"\eta = %$eta" for eta in eta_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(η=η)), Tsit5(), saveat=0.5) for η in η_vals];
+sols = [solve(ODEProblem(F, x_0, tspan, p_gen(eta=eta)), Tsit5(), saveat=0.5) for eta in eta_vals];
 ```
 
 Next, plot the $R_0$ over time:
@@ -468,21 +469,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₀_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)
+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)
 
 
 # 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₀_L, 0.0, 0.0] # start in lockdown
+x_0 = [s_0, e_0, i_0, 0.0, R_0_L, 0.0, 0.0] # start in lockdown
 
-# create two problems, with rapid movement of R₀(t) towards R̄₀(t)
+# create two problems, with rapid movement of R_0(t) towards R_hat_0(t)
 prob_early = ODEProblem(F, x_0, tspan, p_early)
 prob_late = ODEProblem(F, x_0, tspan, p_late)
 ```
@@ -501,7 +502,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.δ * p.γ * sol[3,:]
+flow_deaths(sol, p) = p.N * p.delta * p.gamma * 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")