From 0dbf48d8faca768798e7d8604700ecc6d1fb42df Mon Sep 17 00:00:00 2001 From: annabellasd Date: Mon, 25 Sep 2023 20:20:43 -0800 Subject: [PATCH 1/9] updated lectures --- lectures/continuous_time/covid_sde.md | 121 +++++++++++++------------ lectures/continuous_time/seir_model.md | 81 +++++++++-------- 2 files changed, 103 insertions(+), 99 deletions(-) diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index c0576465..b13b6b85 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/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, trajectories = 100, saveat = 1.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$", + labels = L"Middle 95% Quantile, $\eta = %$eta_1$", layout = (2, 1), size = (900, 900), fillalpha = 0.5) -plot!(summ_2, idxs = (4,5), + +plot!(summ_2, idxs = (4, 5), legend = [:topleft :bottomright], - labels = L"Middle 95% Quantile, $\eta = %$η_2$", size = (900, 900), fillalpha = 0.5) + 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 diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index c36fdb06..14c98d0c 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; γ = 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_0_bar, eta, delta) = p + + return [-gamma * R_0 * s * i; # ds/dt + gamma * R_0 * s * i - sigma * e; # de/dt + sigma * e - γ*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_0_bar = (t, p) -> p.R_0_n) = (T, gamma, sigma, eta, R_0_n, delta, N, R_0_bar) ``` -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,11 +469,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₀_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 @@ -480,9 +481,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₀_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") From ab24bd648465c306e27b99b6e59f7c4099658f4d Mon Sep 17 00:00:00 2001 From: annabellasd Date: Mon, 25 Sep 2023 20:27:31 -0800 Subject: [PATCH 2/9] Update seir_model.md --- lectures/continuous_time/seir_model.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index 14c98d0c..61c5afe8 100644 --- a/lectures/continuous_time/seir_model.md +++ b/lectures/continuous_time/seir_model.md @@ -294,7 +294,7 @@ function F(x, p, t) return [-gamma * R_0 * s * i; # ds/dt gamma * R_0 * s * i - sigma * e; # de/dt - sigma * e - γ*i; # di/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 From f0317ce942254ea90ef9d4ed83d10e8ff117a594 Mon Sep 17 00:00:00 2001 From: annabellasd Date: Thu, 28 Sep 2023 12:50:20 -0800 Subject: [PATCH 3/9] update --- lectures/continuous_time/covid_sde.md | 2 +- lectures/continuous_time/seir_model.md | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index b13b6b85..d1435efa 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -371,7 +371,7 @@ 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, trajectories = 100, saveat = 1.0, x_0, T = 120.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) diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index 61c5afe8..2d17cb2f 100644 --- a/lectures/continuous_time/seir_model.md +++ b/lectures/continuous_time/seir_model.md @@ -290,7 +290,7 @@ 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_0_bar, eta, delta) = p + (;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 @@ -313,9 +313,9 @@ The only confusing part of the notation is the `R_0(t, p)` which evaluates the ` 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, +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_0_bar = (t, p) -> p.R_0_n) = (T, gamma, sigma, eta, R_0_n, delta, N, R_0_bar) + 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_hat_0` function From 06035485c4ff0221ceed11b638e3420ab59ba43b Mon Sep 17 00:00:00 2001 From: annabellasd Date: Thu, 28 Sep 2023 18:04:25 -0800 Subject: [PATCH 4/9] Update covid_sde.md --- lectures/continuous_time/covid_sde.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index d1435efa..0d57d381 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -554,8 +554,8 @@ p_re_gen(; T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20, 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_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_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 From 4b3477ec4c3424c396f98c4571c4460038f75e46 Mon Sep 17 00:00:00 2001 From: annabellasd Date: Wed, 4 Oct 2023 20:36:52 -0800 Subject: [PATCH 5/9] Update covid_sde.md --- lectures/continuous_time/covid_sde.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index 0d57d381..ff1bee0d 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -387,12 +387,12 @@ 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], + 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], + legend = :topleft, labels = L"Middle 95% Quantile, $\eta = %$eta_2$", size = (900, 900), fillalpha = 0.5) ``` From 088e6a29ac5212a71dbd4c8d46319737b3dede3a Mon Sep 17 00:00:00 2001 From: annabellasd Date: Thu, 5 Oct 2023 20:54:04 -0800 Subject: [PATCH 6/9] spacing updates --- lectures/continuous_time/covid_sde.md | 22 +++++++++++----------- lectures/continuous_time/seir_model.md | 14 +++++++------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index ff1bee0d..15d148d1 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -231,20 +231,20 @@ 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 + (;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 + (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 + theta * (delta_bar - delta); # ddelta/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 + (;gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p return [0; 0; 0; 0; sigma * sqrt(R_0); xi * sqrt(delta * (1 - delta))] end @@ -253,9 +253,9 @@ 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, +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) + 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 @@ -549,13 +549,13 @@ function F_reinfect(x, p, t) ] end -p_re_gen(; T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20, +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) + 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_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_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 diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index 2d17cb2f..06606dd7 100644 --- a/lectures/continuous_time/seir_model.md +++ b/lectures/continuous_time/seir_model.md @@ -293,8 +293,8 @@ function F(x, p, t) (;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 * 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 @@ -306,7 +306,7 @@ 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_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_bar_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 @@ -318,7 +318,7 @@ p_gen(;T = 550.0, gamma = 1.0 / 18, sigma = 1 / 5.2, eta = 1.0 / 20, 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_hat_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_bar_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$ @@ -472,8 +472,8 @@ and 75,000 agents already exposed to the virus and thus soon to be contagious. 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) +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 @@ -483,7 +483,7 @@ 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 -# create two problems, with rapid movement of R_0(t) towards R_hat_0(t) +# create two problems, with rapid movement of R_0(t) towards R_bar_0(t) prob_early = ODEProblem(F, x_0, tspan, p_early) prob_late = ODEProblem(F, x_0, tspan, p_late) ``` From befd388cc847d1d11caa0b329a9375f03e07c378 Mon Sep 17 00:00:00 2001 From: annabellasd Date: Tue, 10 Oct 2023 15:34:07 -0800 Subject: [PATCH 7/9] updates --- lectures/multi_agent_models/aiyagari.md | 50 +++--- lectures/multi_agent_models/arellano.md | 69 ++++----- lectures/multi_agent_models/harrison_kreps.md | 22 +-- lectures/multi_agent_models/lake_model.md | 146 +++++++++--------- lectures/multi_agent_models/lucas_model.md | 42 ++--- lectures/multi_agent_models/markov_asset.md | 92 ++++++----- lectures/multi_agent_models/markov_perf.md | 40 ++--- lectures/multi_agent_models/matsuyama.md | 104 ++++++------- .../rational_expectations.md | 32 ++-- lectures/multi_agent_models/schelling.md | 6 +- .../multi_agent_models/uncertainty_traps.md | 89 +++++------ 11 files changed, 343 insertions(+), 349 deletions(-) diff --git a/lectures/multi_agent_models/aiyagari.md b/lectures/multi_agent_models/aiyagari.md index 11ae13bb..de0f8a5e 100644 --- a/lectures/multi_agent_models/aiyagari.md +++ b/lectures/multi_agent_models/aiyagari.md @@ -210,23 +210,23 @@ using LaTeXStrings, Parameters, Plots, QuantEcon ``` ```{code-cell} julia -Household = @with_kw (r = 0.01, - w = 1.0, - σ = 1.0, - β = 0.96, - z_chain = MarkovChain([0.9 0.1; 0.1 0.9], [0.1; 1.0]), - a_min = 1e-10, - a_max = 18.0, - a_size = 200, - a_vals = range(a_min, a_max, length = a_size), - z_size = length(z_chain.state_values), - n = a_size * z_size, - s_vals = gridmake(a_vals, z_chain.state_values), - s_i_vals = gridmake(1:a_size, 1:z_size), - u = σ == 1 ? x -> log(x) : x -> (x^(1 - σ) - 1) / (1 - σ), - R = setup_R!(fill(-Inf, n, a_size), a_vals, s_vals, r, w, u), - # -Inf is the utility of dying (0 consumption) - Q = setup_Q!(zeros(n, a_size, n), s_i_vals, z_chain)) +Household(;r = 0.01, + w = 1.0, + sigma = 1.0, + beta = 0.96, + z_chain = MarkovChain([0.9 0.1; 0.1 0.9], [0.1; 1.0]), + a_min = 1e-10, + a_max = 18.0, + a_size = 200, + a_vals = range(a_min, a_max, length = a_size), + z_size = length(z_chain.state_values), + n = a_size * z_size, + s_vals = gridmake(a_vals, z_chain.state_values), + s_i_vals = gridmake(1:a_size, 1:z_size), + u = sigma == 1 ? x -> log(x) : x -> (x^(1 - sigma) - 1) / (1 - sigma), + R = setup_R!(fill(-Inf, n, a_size), a_vals, s_vals, r, w, u), + # -Inf is the utility of dying (0 consumption) + Q = setup_Q!(zeros(n, a_size, n), s_i_vals, z_chain)) = (; r, w, sigma, beta, z_chain, a_min, a_max, a_size, a_vals, z_size, n, s_vals, s_i_vals, u, R, Q) function setup_Q!(Q, s_i_vals, z_chain) for next_s_i in 1:size(Q, 3) @@ -274,7 +274,7 @@ Random.seed!(42); am = Household(a_max = 20.0, r = 0.03, w = 0.956) # Use the instance to build a discrete dynamic program -am_ddp = DiscreteDP(am.R, am.Q, am.β) +am_ddp = DiscreteDP(am.R, am.Q, am.beta) # Solve using policy function iteration results = solve(am_ddp, PFI) @@ -326,16 +326,16 @@ Random.seed!(42); # Firms' parameters const A = 1 const N = 1 -const α = 0.33 -const β = 0.96 -const δ = 0.05 +const alpha = 0.33 +const beta = 0.96 +const delta = 0.05 function r_to_w(r) - return A * (1 - α) * (A * α / (r + δ)) ^ (α / (1 - α)) + return A * (1 - alpha) * (A * alpha / (r + delta)) ^ (alpha / (1 - alpha)) end function rd(K) - return A * α * (N / K) ^ (1 - α) - δ + return A * alpha * (N / K) ^ (1 - alpha) - delta end function prices_to_capital_stock(am, r) @@ -345,7 +345,7 @@ function prices_to_capital_stock(am, r) (;a_vals, s_vals, u) = am setup_R!(am.R, a_vals, s_vals, r, w, u) - aiyagari_ddp = DiscreteDP(am.R, am.Q, am.β) + aiyagari_ddp = DiscreteDP(am.R, am.Q, am.beta) # Compute the optimal policy results = solve(aiyagari_ddp, PFI) @@ -358,7 +358,7 @@ function prices_to_capital_stock(am, r) end # Create an instance of Household -am = Household(β = β, a_max = 20.0) +am = Household(beta = beta, a_max = 20.0) # Create a grid of r values at which to compute demand and supply of capital r_vals = range(0.005, 0.04, length = 20) diff --git a/lectures/multi_agent_models/arellano.md b/lectures/multi_agent_models/arellano.md index 8b0fb7cb..e9e92a40 100644 --- a/lectures/multi_agent_models/arellano.md +++ b/lectures/multi_agent_models/arellano.md @@ -312,25 +312,25 @@ using Test ``` ```{code-cell} julia -function ArellanoEconomy(;β = .953, - γ = 2., +function ArellanoEconomy(;beta = .953, + gamma = 2., r = 0.017, - ρ = 0.945, - η = 0.025, - θ = 0.282, + rho = 0.945, + eta = 0.025, + theta = 0.282, ny = 21, nB = 251) # create grids Bgrid = collect(range(-.4, .4, length = nB)) - mc = tauchen(ny, ρ, η) - Π = mc.p + mc = tauchen(ny, rho, eta) + Pi = mc.p ygrid = exp.(mc.state_values) ydefgrid = min.(.969 * mean(ygrid), ygrid) # define value functions # notice ordered different than Python to take - # advantage of column major layout of Julia) + # advantage of column major layout of Julia vf = zeros(nB, ny) vd = zeros(1, ny) vc = zeros(nB, ny) @@ -338,13 +338,12 @@ function ArellanoEconomy(;β = .953, q = ones(nB, ny) .* (1 / (1 + r)) defprob = zeros(nB, ny) - return (β = β, γ = γ, r = r, ρ = ρ, η = η, θ = θ, ny = ny, - nB = nB, ygrid = ygrid, ydefgrid = ydefgrid, - Bgrid = Bgrid, Π = Π, vf = vf, vd = vd, vc = vc, - policy = policy, q = q, defprob = defprob) + return (;beta, gamma, r, rho, eta, theta, ny, + nB, ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, + policy, q, defprob) end -u(ae, c) = c^(1 - ae.γ) / (1 - ae.γ) +u(ae, c) = c^(1 - ae.gamma) / (1 - ae.gamma) function one_step_update!(ae, EV, @@ -352,8 +351,8 @@ function one_step_update!(ae, EVc) # unpack stuff - @unpack β, γ, r, ρ, η, θ, ny, nB = ae - @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae + (;beta, gamma, r, rho, eta, theta, ny, nB) = ae + (;ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob) = ae zero_ind = searchsortedfirst(Bgrid, 0.) for iy in 1:ny @@ -361,7 +360,7 @@ function one_step_update!(ae, ydef = ae.ydefgrid[iy] # value of being in default with income y - defval = u(ae, ydef) + β * (θ * EVc[zero_ind, iy] + (1-θ) * EVd[1, iy]) + defval = u(ae, ydef) + beta * (theta * EVc[zero_ind, iy] + (1-theta) * EVd[1, iy]) ae.vd[1, iy] = defval for ib in 1:nB @@ -371,7 +370,7 @@ function one_step_update!(ae, pol_ind = 0 for ib_next=1:nB c = max(y - ae.q[ib_next, iy]*Bgrid[ib_next] + B, 1e-14) - m = u(ae, c) + β * EV[ib_next, iy] + m = u(ae, c) + beta * EV[ib_next, iy] if m > current_max current_max = m @@ -390,14 +389,14 @@ end function compute_prices!(ae) # unpack parameters - @unpack β, γ, r, ρ, η, θ, ny, nB = ae + (;beta, gamma, r, rho, eta, theta, ny, nB) = ae # create default values with a matching size vd_compat = repeat(ae.vd, nB) default_states = vd_compat .> ae.vc # update default probabilities and prices - copyto!(ae.defprob, default_states * ae.Π') + copyto!(ae.defprob, default_states * ae.Pi') copyto!(ae.q, (1 .- ae.defprob) / (1 + r)) return end @@ -405,9 +404,9 @@ end function vfi!(ae; tol = 1e-8, maxit = 10000) # unpack stuff - @unpack β, γ, r, ρ, η, θ, ny, nB = ae - @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae - Πt = Π' + (;beta, gamma, r, rho, eta, theta, ny, nB) = ae + (;ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob) = ae + Pit = Pi' # Iteration stuff it = 0 @@ -420,11 +419,11 @@ function vfi!(ae; tol = 1e-8, maxit = 10000) it += 1 # compute expectations for this iterations - # (we need Π' because of order value function dimensions) + # (we need Pi' because of order value function dimensions) copyto!(V_upd, ae.vf) - EV = ae.vf * Πt - EVd = ae.vd * Πt - EVc = ae.vc * Πt + EV = ae.vf * Pit + EVd = ae.vd * Pit + EVc = ae.vc * Pit # update value function one_step_update!(ae, EV, EVd, EVc) @@ -452,7 +451,7 @@ function QuantEcon.simulate(ae, B_init_ind = searchsortedfirst(ae.Bgrid, B_init) # create a QE MarkovChain - mc = MarkovChain(ae.Π) + mc = MarkovChain(ae.Pi) y_sim_indices = simulate(mc, capT + 1; init = y_init_ind) # allocate and fill output @@ -495,8 +494,8 @@ function QuantEcon.simulate(ae, y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]] q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]] - # with probability θ exit default status - default_status[t + 1] = rand() ≥ ae.θ + # with probability theta exit default status + default_status[t + 1] = rand() >= ae.theta end end @@ -591,12 +590,12 @@ using DataFrames, Plots Compute the value function, policy and equilibrium prices ```{code-cell} julia -ae = ArellanoEconomy(β = .953, # time discount rate - γ = 2., # risk aversion +ae = ArellanoEconomy(beta = .953, # time discount rate + gamma = 2., # risk aversion r = 0.017, # international interest rate - ρ = .945, # persistence in output - η = 0.025, # st dev of output shock - θ = 0.282, # prob of regaining access + rho = .945, # persistence in output + eta = 0.025, # st dev of output shock + theta = 0.282, # prob of regaining access ny = 21, # number of points in y grid nB = 251) # number of points in B grid @@ -627,7 +626,7 @@ q_low = zeros(0) q_high = zeros(0) for i in 1:ae.nB b = ae.Bgrid[i] - if -0.35 ≤ b ≤ 0 # to match fig 3 of Arellano + if -0.35 <= b <= 0 # to match fig 3 of Arellano push!(x, b) push!(q_low, ae.q[i, iy_low]) push!(q_high, ae.q[i, iy_high]) diff --git a/lectures/multi_agent_models/harrison_kreps.md b/lectures/multi_agent_models/harrison_kreps.md index 156ac748..d23d9997 100644 --- a/lectures/multi_agent_models/harrison_kreps.md +++ b/lectures/multi_agent_models/harrison_kreps.md @@ -290,12 +290,12 @@ Here's a function that can be used to compute these values using LinearAlgebra function price_single_beliefs(transition, dividend_payoff; - β=.75) + beta=.75) # First compute inverse piece - imbq_inv = inv(I - β * transition) + imbq_inv = inv(I - beta * transition) # Next compute prices - prices = β * ((imbq_inv * transition) * dividend_payoff) + prices = beta * ((imbq_inv * transition) * dividend_payoff) return prices end @@ -416,7 +416,7 @@ Here's code to solve for $\bar p$, $\hat p_a$ and $\hat p_b$ using the iterative ```{code-cell} julia function price_optimistic_beliefs(transitions, dividend_payoff; - β=.75, max_iter=50000, + beta=.75, max_iter=50000, tol=1e-16) # We will guess an initial price vector of [0, 0] @@ -424,11 +424,11 @@ function price_optimistic_beliefs(transitions, p_old = [10.0,10.0] # We know this is a contraction mapping, so we can iterate to conv - for i ∈ 1:max_iter + for i in 1:max_iter p_old = p_new temp = [maximum((q * p_old) + (q * dividend_payoff)) for q in transitions] - p_new = β * temp + p_new = beta * temp # If we succed in converging, break out of for loop if maximum(sqrt, ((p_new - p_old).^2)) < 1e-12 @@ -437,7 +437,7 @@ function price_optimistic_beliefs(transitions, end temp=[minimum((q * p_old) + (q * dividend_payoff)) for q in transitions] - ptwiddle = β * temp + ptwiddle = beta * temp phat_a = [p_new[1], ptwiddle[2]] phat_b = [ptwiddle[1], p_new[2]] @@ -486,17 +486,17 @@ Here's code to solve for $\check p$ using iteration ```{code-cell} julia function price_pessimistic_beliefs(transitions, dividend_payoff; - β=.75, max_iter=50000, + beta=.75, max_iter=50000, tol=1e-16) # We will guess an initial price vector of [0, 0] p_new = [0, 0] p_old = [10.0, 10.0] # We know this is a contraction mapping, so we can iterate to conv - for i ∈ 1:max_iter + for i in 1:max_iter p_old = p_new temp=[minimum((q * p_old) + (q* dividend_payoff)) for q in transitions] - p_new = β * temp + p_new = beta * temp # If we succed in converging, break out of for loop if maximum(sqrt, ((p_new - p_old).^2)) < 1e-12 @@ -594,7 +594,7 @@ heterogeneous beliefs. opt_beliefs = price_optimistic_beliefs([qa, qb], dividendreturn) labels = ["p_optimistic", "p_hat_a", "p_hat_b"] -for (p, label) ∈ zip(opt_beliefs, labels) +for (p, label) in zip(opt_beliefs, labels) println(label) println(repeat("=", 20)) s0, s1 = round.(p, digits = 2) diff --git a/lectures/multi_agent_models/lake_model.md b/lectures/multi_agent_models/lake_model.md index 1e87f825..1866c63e 100644 --- a/lectures/multi_agent_models/lake_model.md +++ b/lectures/multi_agent_models/lake_model.md @@ -213,20 +213,20 @@ using Test ``` ```{code-cell} julia -LakeModel = @with_kw (λ = 0.283, α = 0.013, b = 0.0124, d = 0.00822) +LakeModel(;lambda = 0.283, alpha = 0.013, b = 0.0124, d = 0.00822) = (;lambda, alpha, b, d) function transition_matrices(lm) - (;λ, α, b, d) = lm + (;lambda, alpha, b, d) = lm g = b - d - A = [(1 - λ) * (1 - d) + b (1 - d) * α + b - (1 - d) * λ (1 - d) * (1 - α)] - Â = A ./ (1 + g) - return (A = A, Â = Â) + A = [(1 - lambda) * (1 - d) + b (1 - d) * alpha + b + (1 - d) * lambda (1 - d) * (1 - alpha)] + A_hat = A ./ (1 + g) + return (A, A_hat) end function rate_steady_state(lm) - (;Â) = transition_matrices(lm) - sol = fixedpoint(x -> Â * x, fill(0.5, 2)) + (;A_hat) = transition_matrices(lm) + sol = fixedpoint(x -> A_hat * x, fill(0.5, 2)) converged(sol) || error("Failed to converge in $(result.iterations) iterations") return sol.zero end @@ -243,12 +243,12 @@ function simulate_stock_path(lm, X0, T) end function simulate_rate_path(lm, x0, T) - (;Â) = transition_matrices(lm) + (;A_hat) = transition_matrices(lm) x_path = zeros(eltype(x0), 2, T) x = copy(x0) for t in 1:T x_path[:, t] = x - x = Â * x + x = A_hat * x end return x_path end @@ -258,24 +258,24 @@ Let's observe these matrices for the baseline model ```{code-cell} julia lm = LakeModel() -A, Â = transition_matrices(lm) +A, A_hat = transition_matrices(lm) A ``` ```{code-cell} julia -Â +A_hat ``` And a revised model ```{code-cell} julia -lm = LakeModel(α = 2.0) -A, Â = transition_matrices(lm) +lm = LakeModel(alpha = 2.0) +A, A_hat = transition_matrices(lm) A ``` ```{code-cell} julia -Â +A_hat ``` ```{code-cell} julia @@ -283,7 +283,7 @@ Â tags: [remove-cell] --- @testset begin - @test lm.α ≈ 2.0 + @test lm.alpha ≈ 2.0 @test A[1][1] ≈ 0.7235062600000001 end ``` @@ -343,8 +343,8 @@ This is the case for our default parameters: ```{code-cell} julia lm = LakeModel() -A, Â = transition_matrices(lm) -e, f = eigvals(Â) +A, A_hat = transition_matrices(lm) +e, f = eigvals(A_hat) abs(e), abs(f) ``` @@ -487,9 +487,9 @@ using QuantEcon, Roots, Random lm = LakeModel(d = 0, b = 0) T = 5000 # Simulation length -(;α, λ) = lm -P = [(1 - λ) λ - α (1 - α)] +(;alpha, lambda) = lm +P = [(1 - lambda) lambda + alpha (1 - alpha)] ``` ```{code-cell} julia @@ -498,9 +498,9 @@ mc = MarkovChain(P, [0; 1]) # 0=unemployed, 1=employed xbar = rate_steady_state(lm) s_path = simulate(mc, T; init=2) -s̄_e = cumsum(s_path) ./ (1:T) -s̄_u = 1 .- s̄_e -s_bars = [s̄_u s̄_e] +s_bar_e = cumsum(s_path) ./ (1:T) +s_bar_u = 1 .- s_bar_e +s_bars = [s_bar_u s_bar_e] plt_unemp = plot(title = "Percent of time unemployed", 1:T, s_bars[:,1],color = :blue, lw = 2, alpha = 0.5, label = "", grid = true) @@ -630,19 +630,18 @@ Following {cite}`davis2006flow`, we set $\alpha$, the hazard rate of leaving emp We will make use of (with some tweaks) the code we wrote in the {doc}`McCall model lecture <../dynamic_programming/mccall_model>`, embedded below for convenience. ```{code-cell} julia -function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1e-5, - iter = 2_000) - (;α, β, σ, c, γ, w, E, u) = mcm +function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1e-5,iter = 2_000) + (;alpha, beta, sigma, c, gamma, w, E, u) = mcm # necessary objects - u_w = u.(w, σ) - u_c = u(c, σ) + u_w = u.(w, sigma) + u_c = u(c, sigma) # Bellman operator T. Fixed point is x* s.t. T(x*) = x* function T(x) V = x[1:end-1] U = x[end] - [u_w + β * ((1 - α) * V .+ α * U); u_c + β * (1 - γ) * U + β * γ * E * max.(U, V)] + [u_w + beta * ((1 - alpha) * V .+ alpha * U); u_c + beta * (1 - gamma) * U + beta * gamma * E * max.(U, V)] end # value function iteration @@ -654,13 +653,13 @@ function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1 # compute the reservation wage w_barindex = searchsortedfirst(V .- U, 0.0) if w_barindex >= length(w) # if this is true, you never want to accept - w̄ = Inf + w_bar = Inf else - w̄ = w[w_barindex] # otherwise, return the number + w_bar = w[w_barindex] # otherwise, return the number end # return a NamedTuple, so we can select values by name - return (V = V, U = U, w̄ = w̄) + return (;V, U, w_bar) end ``` @@ -668,17 +667,18 @@ And the McCall object ```{code-cell} julia # a default utility function -u(c, σ) = c > 0 ? (c^(1 - σ) - 1) / (1 - σ) : -10e-6 +u(c, sigma) = c > 0 ? (c^(1 - sigma) - 1) / (1 - sigma) : -10e-6 # model constructor -McCallModel = @with_kw (α = 0.2, - β = 0.98, # discount rate - γ = 0.7, - c = 6.0, # unemployment compensation - σ = 2.0, - u = u, # utility function - w = range(10, 20, length = 60), # wage values - E = Expectation(BetaBinomial(59, 600, 400))) # distribution over wage values +McCallModel(;alpha = 0.2, + beta = 0.98, # discount rate + gamma = 0.7, + c = 6.0, # unemployment compensation + sigma = 2.0, + u = u, # utility function + w = range(10, 20, length = 60), # w= wage values; E= distribution over wage values + E = Expectation(BetaBinomial(59, 600, 400))) = (;alpha, beta, gamma, c, sigma, u, w, E) + ``` Now let's compute and plot welfare, employment, unemployment, and tax revenue as a @@ -686,13 +686,13 @@ function of the unemployment compensation rate ```{code-cell} julia # some global variables that will stay constant -α = 0.013 -α_q = (1 - (1 - α)^3) +alpha = 0.013 +alpha_q = (1 - (1 - alpha)^3) b_param = 0.0124 d_param = 0.00822 -β = 0.98 -γ = 1.0 -σ = 2.0 +beta = 0.98 +gamma = 1.0 +sigma = 2.0 # the default wage distribution: a discretized log normal log_wage_mean, wage_grid_size, max_wage = 20, 200, 170 @@ -707,33 +707,33 @@ w_vec = (w_vec[1:end-1] + w_vec[2:end]) / 2 E = expectation(Categorical(p_vec)) # expectation object -function compute_optimal_quantities(c, τ) - mcm = McCallModel(α = α_q, - β = β, - γ = γ, - c = c - τ, # post-tax compensation - σ = σ, - w = w_vec .- τ, # post-tax wages +function compute_optimal_quantities(c, tau) + mcm = McCallModel(alpha = alpha_q, + beta = beta, + gamma = gamma, + c = c - tau, # post-tax compensation + sigma = sigma, + w = w_vec .- tau, # post-tax wages E = E) # expectation operator - (;V, U, w̄) = solve_mccall_model(mcm) - indicator = wage -> wage > w̄ - λ = γ * E * indicator.(w_vec .- τ) + (;V, U, w_bar) = solve_mccall_model(mcm) + indicator = wage -> wage > w_bar + lambda = gamma * E * indicator.(w_vec .- tau) - return w̄, λ, V, U + return w_bar, lambda, V, U end -function compute_steady_state_quantities(c, τ) - w̄, λ_param, V, U = compute_optimal_quantities(c, τ) +function compute_steady_state_quantities(c, tau) + w_bar, lambda_param, V, U = compute_optimal_quantities(c, tau) # compute steady state employment and unemployment rates - lm = LakeModel(λ = λ_param, α = α_q, b = b_param, d = d_param) + lm = LakeModel(lambda = lambda_param, alpha = alpha_q, b = b_param, d = d_param) x = rate_steady_state(lm) u_rate, e_rate = x # compute steady state welfare - indicator(wage) = wage > w̄ - decisions = indicator.(w_vec .- τ) + indicator(wage) = wage > w_bar + decisions = indicator.(w_vec .- tau) w = (E * (V .* decisions)) / (E * decisions) welfare = e_rate .* w + u_rate .* U @@ -746,8 +746,8 @@ function find_balanced_budget_tax(c) return t - u_rate * c end - τ = find_zero(steady_state_budget, (0.0, 0.9c)) - return τ + tau = find_zero(steady_state_budget, (0.0, 0.9c)) + return tau end # levels of unemployment insurance we wish to study @@ -844,7 +844,7 @@ T = 50 New legislation changes $\lambda$ to $0.2$ ```{code-cell} julia -lm = LakeModel(λ = 0.2) +lm = LakeModel(lambda = 0.2) ``` ```{code-cell} julia @@ -936,16 +936,16 @@ end Here are the other parameters: ```{code-cell} julia -b̂ = 0.003 -T̂ = 20 +b_hat = 0.003 +T_hat = 20 ``` Let's increase $b$ to the new value and simulate for 20 periods ```{code-cell} julia -lm = LakeModel(b=b̂) -X_path1 = simulate_stock_path(lm, x0 * N0, T̂) # simulate stocks -x_path1 = simulate_rate_path(lm, x0, T̂) # simulate rates +lm = LakeModel(b=b_hat) +X_path1 = simulate_stock_path(lm, x0 * N0, T_hat) # simulate stocks +x_path1 = simulate_rate_path(lm, x0, T_hat) # simulate rates ``` Now we reset $b$ to the original value and then, using the state @@ -954,8 +954,8 @@ additional 30 periods ```{code-cell} julia lm = LakeModel(b = 0.0124) -X_path2 = simulate_stock_path(lm, X_path1[:, end-1], T-T̂+1) # simulate stocks -x_path2 = simulate_rate_path(lm, x_path1[:, end-1], T-T̂+1) # simulate rates +X_path2 = simulate_stock_path(lm, X_path1[:, end-1], T-T_hat+1) # simulate stocks +x_path2 = simulate_rate_path(lm, x_path1[:, end-1], T-T_hat+1) # simulate rates ``` Finally we combine these two paths and plot diff --git a/lectures/multi_agent_models/lucas_model.md b/lectures/multi_agent_models/lucas_model.md index b5239a57..2b26c385 100644 --- a/lectures/multi_agent_models/lucas_model.md +++ b/lectures/multi_agent_models/lucas_model.md @@ -393,39 +393,39 @@ using Distributions, Interpolations, LaTeXStrings, Parameters, Plots, Random ```{code-cell} julia # model -function LucasTree(;γ = 2.0, - β = 0.95, - α = 0.9, - σ = 0.1, +function LucasTree(;gamma = 2.0, + beta = 0.95, + alpha = 0.9, + sigma = 0.1, grid_size = 100) - ϕ = LogNormal(0.0, σ) - shocks = rand(ϕ, 500) + phi = LogNormal(0.0, sigma) + shocks = rand(phi, 500) # build a grid with mass around stationary distribution - ssd = σ / sqrt(1 - α^2) + ssd = sigma / sqrt(1 - alpha^2) grid_min, grid_max = exp(-4ssd), exp(4ssd) grid = range(grid_min, grid_max, length = grid_size) - # set h(y) = β * int u'(G(y,z)) G(y,z) ϕ(dz) + # set h(y) = beta * int u'(G(y,z)) G(y,z) phi(dz) h = similar(grid) for (i, y) in enumerate(grid) - h[i] = β * mean((y^α .* shocks).^(1 - γ)) + h[i] = beta * mean((y^alpha .* shocks).^(1 - gamma)) end - return (γ = γ, β = β, α = α, σ = σ, ϕ = ϕ, grid = grid, shocks = shocks, h = h) + return (gamma = gamma, beta = beta, alpha = alpha, sigma = sigma, phi = phi, grid = grid, shocks = shocks, h = h) end # approximate Lucas operator, which returns the updated function Tf on the grid function lucas_operator(lt, f) # unpack input - (;grid, α, β, h) = lt + (;grid, alpha, beta, h) = lt z = lt.shocks - Af = LinearInterpolation(grid, f, extrapolation_bc=Line()) + Af = linear_interpolation(grid, f, extrapolation_bc=Line()) - Tf = [ h[i] + β * mean(Af.(grid[i]^α .* z)) for i in 1:length(grid) ] + Tf = [ h[i] + beta * mean(Af.(grid[i]^alpha .* z)) for i in 1:length(grid) ] return Tf end @@ -434,7 +434,7 @@ function solve_lucas_model(lt; tol = 1e-6, max_iter = 500) - (;grid, γ) = lt + (;grid, gamma) = lt i = 0 f = zero(grid) # Initial guess of f @@ -447,8 +447,8 @@ function solve_lucas_model(lt; i += 1 end - # p(y) = f(y) * y ^ γ - price = f .* grid.^γ + # p(y) = f(y) * y ^ gamma + price = f .* grid.^gamma return price end @@ -459,7 +459,7 @@ An example of usage is given in the docstring and repeated here ```{code-cell} julia Random.seed!(42) # For reproducible results. -tree = LucasTree(γ = 2.0, β = 0.95, α = 0.90, σ = 0.1) +tree = LucasTree(gamma = 2.0, beta = 0.95, alpha = 0.90, sigma = 0.1) price_vals = solve_lucas_model(tree); ``` @@ -519,11 +519,11 @@ Random.seed!(42); ```{code-cell} julia plot() -for β in (.95, 0.98) - tree = LucasTree(;β = β) +for beta in (.95, 0.98) + tree = LucasTree(;beta = beta) grid = tree.grid price_vals = solve_lucas_model(tree) - plot!(grid, price_vals, lw = 2, label = L"\beta = %$β") + plot!(grid, price_vals, lw = 2, label = L"\beta = %$beta") end plot!(xlabel = L"y", ylabel = "price", legend = :topleft) @@ -535,7 +535,7 @@ tags: [remove-cell] --- @testset begin # For the 0.98, since the other one is overwritten. Random.seed!(42) - price_vals = solve_lucas_model(LucasTree(β = 0.98)) + price_vals = solve_lucas_model(LucasTree(beta = 0.98)) #test price_vals[20] ≈ 35.00073581199659 #test price_vals[57] ≈ 124.32987344509688 end diff --git a/lectures/multi_agent_models/markov_asset.md b/lectures/multi_agent_models/markov_asset.md index 27bc9878..821c702c 100644 --- a/lectures/multi_agent_models/markov_asset.md +++ b/lectures/multi_agent_models/markov_asset.md @@ -404,12 +404,12 @@ Here's the code, including a test of the spectral radius condition ```{code-cell} julia n = 25 # size of state space -β = 0.9 +beta = 0.9 mc = tauchen(n, 0.96, 0.02) K = mc.p .* exp.(mc.state_values)' -v = (I - β * K) \ (β * K * ones(n, 1)) +v = (I - beta * K) \ (beta * K * ones(n, 1)) plot(mc.state_values, v, @@ -541,39 +541,35 @@ the AssetPriceModel objects ```{code-cell} julia # A default Markov chain for the state process -ρ = 0.9 -σ = 0.02 +rho = 0.9 +sigma = 0.02 n = 25 -default_mc = tauchen(n, ρ, σ) +default_mc = tauchen(n, rho, sigma) -AssetPriceModel = @with_kw (β = 0.96, - γ = 2.0, - mc = default_mc, - n = size(mc.p)[1], - g = exp) +AssetPriceModel(;beta = 0.96, gamma = 2.0, mc = default_mc, n = size(mc.p)[1], g = exp) = (;beta, gamma, mc, n, g) # test stability of matrix Q function test_stability(ap, Q) sr = maximum(abs, eigvals(Q)) - if sr ≥ 1 / ap.β + if sr >= 1 / ap.beta msg = "Spectral radius condition failed with radius = $sr" throw(ArgumentError(msg)) end end # price/dividend ratio of the Lucas tree -function tree_price(ap; γ = ap.γ) +function tree_price(ap; gamma = ap.gamma) # Simplify names, set up matrices - (;β, mc) = ap + (;beta, mc) = ap P, y = mc.p, mc.state_values y = reshape(y, 1, ap.n) - J = P .* ap.g.(y).^(1 - γ) + J = P .* ap.g.(y).^(1 - gamma) # Make sure that a unique solution exists test_stability(ap, J) # Compute v - v = (I - β * J) \ sum(β * J, dims = 2) + v = (I - beta * J) \ sum(beta * J, dims = 2) return v end @@ -583,16 +579,16 @@ Here's a plot of $v$ as a function of the state for several values of $\gamma$, with a positively correlated Markov process and $g(x) = \exp(x)$ ```{code-cell} julia -γs = [1.2, 1.4, 1.6, 1.8, 2.0] +gammas = [1.2, 1.4, 1.6, 1.8, 2.0] ap = AssetPriceModel() states = ap.mc.state_values lines = [] labels = [] -for γ in γs - v = tree_price(ap, γ = γ) - label = L"\gamma = %$γ" +for gamma in gammas + v = tree_price(ap, gamma = gamma) + label = L"\gamma = %$gamma" push!(labels, label) push!(lines, v) end @@ -690,18 +686,18 @@ p = (I - \beta M)^{-1} \beta M \zeta {\mathbb 1} The above is implemented in the function consol_price ```{code-cell} julia -function consol_price(ap, ζ) +function consol_price(ap, zeta) # Simplify names, set up matrices - (;β, γ, mc, g, n) = ap + (;beta, gamma, mc, g, n) = ap P, y = mc.p, mc.state_values y = reshape(y, 1, n) - M = P .* g.(y).^(-γ) + M = P .* g.(y).^(-gamma) # Make sure that a unique solution exists test_stability(ap, M) # Compute price - return (I - β * M) \ sum(β * ζ * M, dims = 2) + return (I - beta * M) \ sum(beta * zeta * M, dims = 2) end ``` @@ -777,24 +773,24 @@ We can find the solution with the following function call_option ```{code-cell} julia # price of perpetual call on consol bond -function call_option(ap, ζ, p_s, ϵ = 1e-7) +function call_option(ap, zeta, p_s, epsilon = 1e-7) # Simplify names, set up matrices - (;β, γ, mc, g, n) = ap + (;beta, gamma, mc, g, n) = ap P, y = mc.p, mc.state_values y = reshape(y, 1, n) - M = P .* g.(y).^(-γ) + M = P .* g.(y).^(-gamma) # Make sure that a unique console price exists test_stability(ap, M) # Compute option price - p = consol_price(ap, ζ) + p = consol_price(ap, zeta) w = zeros(ap.n, 1) - error = ϵ + 1 - while (error > ϵ) + error = epsilon + 1 + while (error > epsilon) # Maximize across columns - w_new = max.(β * M * w, p .- p_s) + w_new = max.(beta * M * w, p .- p_s) # Find maximal difference of each component and update error = maximum(abs, w - w_new) w = w_new @@ -807,13 +803,13 @@ end Here's a plot of $w$ compared to the consol price when $P_S = 40$ ```{code-cell} julia -ap = AssetPriceModel(β=0.9) -ζ = 1.0 +ap = AssetPriceModel(beta=0.9) +zeta = 1.0 strike_price = 40.0 x = ap.mc.state_values -p = consol_price(ap, ζ) -w = call_option(ap, ζ, strike_price) +p = consol_price(ap, zeta) +w = call_option(ap, zeta, strike_price) plot(x, p, color = "blue", lw = 2, xlabel = "state", label = "consol price") plot!(x, w, color = "green", lw = 2, label = "value of call option") @@ -899,9 +895,9 @@ Consider the following primitives n = 5 P = fill(0.0125, n, n) + (0.95 - 0.0125)I s = [1.05, 1.025, 1.0, 0.975, 0.95] -γ = 2.0 -β = 0.94 -ζ = 1.0 +gamma = 2.0 +beta = 0.94 +zeta = 1.0 ``` Let $g$ be defined by $g(x) = x$ (that is, $g$ is the identity map). @@ -965,16 +961,16 @@ P = fill(0.0125, n, n) + (0.95 - 0.0125)I s = [0.95, 0.975, 1.0, 1.025, 1.05] # state values mc = MarkovChain(P, s) -γ = 2.0 -β = 0.94 -ζ = 1.0 +gamma = 2.0 +beta = 0.94 +zeta = 1.0 p_s = 150.0 ``` Next we'll create an instance of AssetPriceModel to feed into the functions. ```{code-cell} julia -ap = AssetPriceModel(β = β, mc = mc, γ = γ, g = x -> x) +ap = AssetPriceModel(beta = beta, mc = mc, gamma = gamma, g = x -> x) ``` ```{code-cell} julia @@ -997,7 +993,7 @@ println("Consol Bond Prices: $(v_consol)\n") ``` ```{code-cell} julia -w = call_option(ap, ζ, p_s) +w = call_option(ap, zeta, p_s) ``` ```{code-cell} julia @@ -1015,23 +1011,23 @@ end Here's a suitable function: ```{code-cell} julia -function finite_horizon_call_option(ap, ζ, p_s, k) +function finite_horizon_call_option(ap, zeta, p_s, k) # Simplify names, set up matrices - (;β, γ, mc) = ap + (;beta, gamma, mc) = ap P, y = mc.p, mc.state_values y = y' - M = P .* ap.g.(y).^(- γ) + M = P .* ap.g.(y).^(- gamma) # Make sure that a unique console price exists test_stability(ap, M) # Compute option price - p = consol_price(ap, ζ) + p = consol_price(ap, zeta) w = zeros(ap.n, 1) for i in 1:k # Maximize across columns - w = max.(β * M * w, p .- p_s) + w = max.(beta * M * w, p .- p_s) end return w @@ -1052,7 +1048,7 @@ end lines = [] labels = [] for k in [5, 25] - w = finite_horizon_call_option(ap, ζ, p_s, k) + w = finite_horizon_call_option(ap, zeta, p_s, k) push!(lines, w) push!(labels, L"k = %$k") end diff --git a/lectures/multi_agent_models/markov_perf.md b/lectures/multi_agent_models/markov_perf.md index a7b61d10..06fcfcf9 100644 --- a/lectures/multi_agent_models/markov_perf.md +++ b/lectures/multi_agent_models/markov_perf.md @@ -432,8 +432,8 @@ using QuantEcon, LinearAlgebra # parameters a0 = 10.0 a1 = 2.0 -β = 0.96 -γ = 12.0 +beta = 0.96 +gamma = 12.0 # in LQ form A = I + zeros(3, 3) @@ -448,12 +448,12 @@ R2 = [ 0.0 0.0 -a0 / 2.0; 0.0 0.0 a1 / 2.0; -a0 / 2.0 a1 / 2.0 a1] -Q1 = Q2 = γ +Q1 = Q2 = gamma S1 = S2 = W1 = W2 = M1 = M2 = 0.0 # solve using QE's nnash function F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2, - beta=β) + beta=beta) # display policies println("Computed policies for firm 1 and firm 2:") @@ -480,8 +480,8 @@ In particular, let's take F2 as computed above, plug it into {eq}`eq_mpe_p1p` an We hope that the resulting policy will agree with F1 as computed above ```{code-cell} julia -Λ1 = A - (B2 * F2) -lq1 = QuantEcon.LQ(Q1, R1, Λ1, B1, bet=β) +Lamda1 = A - (B2 * F2) +lq1 = QuantEcon.LQ(Q1, R1, Lamda1, B1, bet=beta) P1_ih, F1_ih, d = stationary_values(lq1) F1_ih ``` @@ -493,7 +493,7 @@ tags: [remove-cell] @testset begin @test P1_ih[2, 2] ≈ 5.441368459897164 @test d ≈ 0.0 - @test Λ1[1, 1] ≈ 1.0 && Λ1[3, 2] ≈ -0.07584666305807419 + @test Lamda1[1, 1] ≈ 1.0 && Lamda1[3, 2] ≈ -0.07584666305807419 @test F1_ih ≈ [-0.6684661291052371 0.29512481789806305 0.07584666292394007] @test isapprox(F1, F1_ih, atol=1e-7) # Make sure the test below comes up true. end @@ -651,7 +651,7 @@ The exercise is to calculate these matrices and compute the following figures. The first figure shows the dynamics of inventories for each firm when the parameters are ```{code-cell} julia -δ = 0.02 +delta = 0.02 D = [ -1 0.5; 0.5 -1] b = [25, 25] @@ -683,8 +683,8 @@ First let's compute the duopoly MPE under the stated parameters # parameters a0 = 10.0 a1 = 2.0 -β = 0.96 -γ = 12.0 +beta = 0.96 +gamma = 12.0 # in LQ form A = I + zeros(3, 3) @@ -699,12 +699,12 @@ R2 = [ 0.0 0.0 -a0 / 2.0; 0.0 0.0 a1 / 2.0; -a0 / 2.0 a1 / 2.0 a1] -Q1 = Q2 = γ +Q1 = Q2 = gamma S1 = S2 = W1 = W2 = M1 = M2 = 0.0 # solve using QE's nnash function F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2, - beta=β) + beta=beta) ``` ```{code-cell} julia @@ -783,9 +783,9 @@ resulting dynamics of $\{q_t\}$, starting at $q_0 = 2.0$. tags: [hide-output] --- R = a1 -Q = γ +Q = gamma A = B = 1 -lq_alt = QuantEcon.LQ(Q, R, A, B, bet=β) +lq_alt = QuantEcon.LQ(Q, R, A, B, bet=beta) P, F, d = stationary_values(lq_alt) q̄ = a0 / (2.0 * a1) qm = zeros(n) @@ -830,13 +830,13 @@ plot(plt_q, plt_p, layout=(2,1), size=(700,600)) We treat the case $\delta = 0.02$ ```{code-cell} julia -δ = 0.02 +delta = 0.02 D = [-1 0.5; 0.5 -1] b = [25, 25] c1 = c2 = [1, -2, 1] e1 = e2 = [10, 10, 3] -δ_1 = 1-δ +delta_1 = 1-delta ``` Recalling that the control and state are @@ -860,14 +860,14 @@ we set up the matrices as follows: ```{code-cell} julia # create matrices needed to compute the Nash feedback equilibrium -A = [δ_1 0 -δ_1 * b[1]; - 0 δ_1 -δ_1 * b[2]; +A = [delta_1 0 -delta_1 * b[1]; + 0 delta_1 -delta_1 * b[2]; 0 0 1] -B1 = δ_1 * [1 -D[1, 1]; +B1 = delta_1 * [1 -D[1, 1]; 0 -D[2, 1]; 0 0] -B2 = δ_1 * [0 -D[1, 2]; +B2 = delta_1 * [0 -D[1, 2]; 1 -D[2, 2]; 0 0] diff --git a/lectures/multi_agent_models/matsuyama.md b/lectures/multi_agent_models/matsuyama.md index 3153335a..88ee3b50 100644 --- a/lectures/multi_agent_models/matsuyama.md +++ b/lectures/multi_agent_models/matsuyama.md @@ -342,7 +342,7 @@ using Plots, Parameters ``` ```{code-cell} julia -function h_j(j, nk, s1, s2, θ, δ, ρ) +function h_j(j, nk, s1, s2, theta, delta, rho) # Find out who's h we are evaluating if j == 1 sj = s1 @@ -354,8 +354,8 @@ function h_j(j, nk, s1, s2, θ, δ, ρ) # Coefficients on the quadratic a x^2 + b x + c = 0 a = 1.0 - b = ((ρ + 1 / ρ) * nk - sj - sk) - c = (nk * nk - (sj * nk) / ρ - sk * ρ * nk) + b = ((rho + 1 / rho) * nk - sj - sk) + c = (nk * nk - (sj * nk) / rho - sk * rho * nk) # Positive solution of quadratic form root = (-b + sqrt(b * b - 4 * a * c)) / (2 * a) @@ -363,41 +363,41 @@ function h_j(j, nk, s1, s2, θ, δ, ρ) return root end -DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≤ s1_ρ) && (n2 ≤ s2_ρ) +DLL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 <= s1_rho) && (n2 <= s2_rho) -DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≥ h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 ≥ h_j(2, n1, s1, s2, θ, δ, ρ)) +DHH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 >= h_j(1, n2, s1, s2, theta, delta, rho)) && (n2 >= h_j(2, n1, s1, s2, theta, delta, rho)) -DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≥ s1_ρ) && (n2 ≤ h_j(2, n1, s1, s2, θ, δ, ρ)) +DHL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 >= s1_rho) && (n2 <= h_j(2, n1, s1, s2, theta, delta, rho)) -DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≤ h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 ≥ s2_ρ) +DLH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 <= h_j(1, n2, s1, s2, theta, delta, rho)) && (n2 >= s2_rho) -function one_step(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) +function one_step(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) # Depending on where we are, evaluate the right branch - if DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * (θ * s1_ρ + (1 - θ) * n1) - n2_tp1 = δ * (θ * s2_ρ + (1 - θ) * n2) - elseif DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * n1 - n2_tp1 = δ * n2 - elseif DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * n1 - n2_tp1 = δ * (θ * h_j(2, n1, s1, s2, θ, δ, ρ) + (1 - θ) * n2) - elseif DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * (θ * h_j(1, n2, s1, s2, θ, δ, ρ) + (1 - θ) * n1) - n2_tp1 = δ * n2 + if DLL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * (theta * s1_rho + (1 - theta) * n1) + n2_tp1 = delta * (theta * s2_rho + (1 - theta) * n2) + elseif DHH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * n1 + n2_tp1 = delta * n2 + elseif DHL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * n1 + n2_tp1 = delta * (theta * h_j(2, n1, s1, s2, theta, delta, rho) + (1 - theta) * n2) + elseif DLH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * (theta * h_j(1, n2, s1, s2, theta, delta, rho) + (1 - theta) * n1) + n2_tp1 = delta * n2 end return n1_tp1, n2_tp1 end -new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - one_step(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) +new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + one_step(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) -function pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, +function pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers) # Initialize the status of synchronization @@ -407,11 +407,11 @@ function pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, nsync = 0 - while (~synchronized) && (iters < maxiter) + while (!synchronized) && (iters < maxiter) # Increment the number of iterations and get next values iters += 1 - n1_t, n2_t = new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) + n1_t, n2_t = new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) # Check whether same in this period if abs(n1_t - n2_t) < 1e-8 @@ -432,7 +432,7 @@ function pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, return synchronized, pers_2_sync end -function create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, +function create_attraction_basis(s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers, npts) # Create unit range with npts synchronized, pers_2_sync = false, 0 @@ -443,9 +443,7 @@ function create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, # Iterate over initial conditions for (i, n1_0) in enumerate(unit_range) for (j, n2_0) in enumerate(unit_range) - synchronized, pers_2_sync = pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, - s1, s2, θ, δ, ρ, - maxiter, npers) + synchronized, pers_2_sync = pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers) time_2_sync[i, j] = pers_2_sync end end @@ -453,18 +451,18 @@ function create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, end # model -function MSGSync(s1 = 0.5, θ = 2.5, δ = 0.7, ρ = 0.2) +function MSGSync(s1 = 0.5, theta = 2.5, delta = 0.7, rho = 0.2) # Store other cutoffs and parameters we use s2 = 1 - s1 - s1_ρ = min((s1 - ρ * s2) / (1 - ρ), 1) - s2_ρ = 1 - s1_ρ + s1_rho = min((s1 - rho * s2) / (1 - rho), 1) + s2_rho = 1 - s1_rho - return (s1 = s1, s2 = s2, s1_ρ = s1_ρ, s2_ρ = s2_ρ, θ = θ, δ = δ, ρ = ρ) + return (;s1, s2, s1_rho, s2_rho, theta, delta, rho) end function simulate_n(model, n1_0, n2_0, T) # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model + (;s1, s2, theta, delta, rho, s1_rho, s2_rho) = model # Allocate space n1 = zeros(T) @@ -474,7 +472,7 @@ function simulate_n(model, n1_0, n2_0, T) for t in 1:T # Get next values n1[t], n2[t] = n1_0, n2_0 - n1_0, n2_0 = new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) + n1_0, n2_0 = new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) end return n1, n2 @@ -483,10 +481,10 @@ end function pers_till_sync(model, n1_0, n2_0, maxiter = 500, npers = 3) # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model + (;s1, s2, theta, delta, rho, s1_rho, s2_rho) = model - return pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, - θ, δ, ρ, maxiter, npers) + return pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, + theta, delta, rho, maxiter, npers) end function create_attraction_basis(model; @@ -494,10 +492,10 @@ function create_attraction_basis(model; npers = 3, npts = 50) # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model + (;s1, s2, theta, delta, rho, s1_rho, s2_rho) = model - ab = create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, - ρ, maxiter, npers, npts) + ab = create_attraction_basis(s1_rho, s2_rho, s1, s2, theta, delta, + rho, maxiter, npers, npts) return ab end ``` @@ -513,8 +511,8 @@ The time series share parameters but differ in their initial condition. Here's the function ```{code-cell} julia -function plot_timeseries(n1_0, n2_0, s1 = 0.5, θ = 2.5, δ = 0.7, ρ = 0.2) - model = MSGSync(s1, θ, δ, ρ) +function plot_timeseries(n1_0, n2_0, s1 = 0.5, theta = 2.5, delta = 0.7, rho = 0.2) + model = MSGSync(s1, theta, delta, rho) n1, n2 = simulate_n(model, n1_0, n2_0, 25) return [n1 n2] end @@ -569,10 +567,10 @@ Replicate the figure {ref}`shown above ` by coloring initial conditions ### Exercise 1 ```{code-cell} julia -function plot_attraction_basis(s1 = 0.5, θ = 2.5, δ = 0.7, ρ = 0.2; npts = 250) +function plot_attraction_basis(s1 = 0.5, theta = 2.5, delta = 0.7, rho = 0.2; npts = 250) # Create attraction basis unitrange = range(0, 1, length = npts) - model = MSGSync(s1, θ, δ, ρ) + model = MSGSync(s1, theta, delta, rho) ab = create_attraction_basis(model,npts=npts) plt = Plots.heatmap(ab, legend = false) end @@ -593,13 +591,13 @@ plot(plots..., size = (1000, 1000)) tags: [remove-cell] --- function plot_attraction_basis(s1 = 0.5, - θ = 2.5, - δ = 0.7, - ρ = 0.2; + theta = 2.5, + delta = 0.7, + rho = 0.2; npts = 250) # Create attraction basis unitrange = range(0, 1, length = npts) - model = MSGSync(s1, θ, δ, ρ) + model = MSGSync(s1, theta, delta, rho) ab = create_attraction_basis(model, npts = npts) end diff --git a/lectures/multi_agent_models/rational_expectations.md b/lectures/multi_agent_models/rational_expectations.md index 4c3d9fe1..2215b90e 100644 --- a/lectures/multi_agent_models/rational_expectations.md +++ b/lectures/multi_agent_models/rational_expectations.md @@ -662,16 +662,16 @@ Here's our solution # model parameters a0 = 100 a1 = 0.05 -β = 0.95 -γ = 10.0 +beta = 0.95 +gamma = 10.0 # beliefs -κ0 = 95.5 -κ1 = 0.95 +kappa0 = 95.5 +kappa1 = 0.95 # formulate the LQ problem A = [1 0 0 - 0 κ1 κ0 + 0 kappa1 kappa0 0 0 1] B = [1.0, 0.0, 0.0] @@ -680,10 +680,10 @@ R = [ 0 a1/2 -a0/2 a1/2 0 0 -a0/2 0 0] -Q = 0.5 * γ +Q = 0.5 * gamma # solve for the optimal policy -lq = QuantEcon.LQ(Q, R, A, B; bet = β) +lq = QuantEcon.LQ(Q, R, A, B; bet = beta) P, F, d = stationary_values(lq) hh = h0, h1, h2 = -F[3], 1 - F[1], -F[2] @@ -751,7 +751,7 @@ for (k0, k1) in candidates A = [1 0 0 0 k1 k0 0 0 1] - lq = QuantEcon.LQ(Q, R, A, B; bet=β) + lq = QuantEcon.LQ(Q, R, A, B; bet=beta) P, F, d = stationary_values(lq) hh = h0, h1, h2 = -F[3], 1 - F[1], -F[2] @@ -833,15 +833,15 @@ B = [1.0, 0.0] R = [ a1 / 2.0 -a0 / 2.0 -a0 / 2.0 0.0] -Q = γ / 2.0 +Q = gamma / 2.0 # solve for the optimal policy -lq = QuantEcon.LQ(Q, R, A, B; bet=β) +lq = QuantEcon.LQ(Q, R, A, B; bet=beta) P, F, d = stationary_values(lq) # print the results -κ0, κ1 = -F[2], 1 - F[1] -println("κ0=$κ0\tκ1=$κ1") +kappa0, kappa1 = -F[2], 1 - F[1] +println("kappa0=$kappa0\tkappa1=$kappa1") ``` ```{code-cell} julia @@ -849,8 +849,8 @@ println("κ0=$κ0\tκ1=$κ1") tags: [remove-cell] --- @testset begin - @test κ0 ≈ 95.081845248600 rtol = 4 - @test κ1 ≈ 0.952459076301 rtol = 4 + @test kappa0 ≈ 95.081845248600 rtol = 4 + @test kappa1 ≈ 0.952459076301 rtol = 4 end ``` @@ -879,10 +879,10 @@ B = [1.0, 0.0] R = [ a1 -a0 / 2.0 -a0 / 2.0 0.0] -Q = γ / 2.0 +Q = gamma / 2.0 # solve for the optimal policy -lq = QuantEcon.LQ(Q, R, A, B; bet=β) +lq = QuantEcon.LQ(Q, R, A, B; bet=beta) P, F, d = stationary_values(lq) # print the results diff --git a/lectures/multi_agent_models/schelling.md b/lectures/multi_agent_models/schelling.md index 868aed51..4bb0f874 100644 --- a/lectures/multi_agent_models/schelling.md +++ b/lectures/multi_agent_models/schelling.md @@ -8,7 +8,7 @@ kernelspec: language: julia name: julia-1.9 --- - +##AT THIS ONE BUT LUCAS (schelling)= ```{raw} html
@@ -160,7 +160,7 @@ Random.seed!(42); ``` ```{code-cell} julia -Agent = @with_kw (kind, location = rand(2)) +Agent(;kind, location = rand(2)) = (;kind, location) draw_location!(a) = a.location .= rand(2) @@ -179,7 +179,7 @@ function is_happy(a) # partialsortperm(get_distance.(Ref(a), agents), # 1:neighborhood_size)) - return share ≥ preference + return share >= preference end function update!(a) diff --git a/lectures/multi_agent_models/uncertainty_traps.md b/lectures/multi_agent_models/uncertainty_traps.md index a637b118..477ef665 100644 --- a/lectures/multi_agent_models/uncertainty_traps.md +++ b/lectures/multi_agent_models/uncertainty_traps.md @@ -236,17 +236,18 @@ using DataFrames, LaTeXStrings, Parameters, Plots ``` ```{code-cell} julia -UncertaintyTrapEcon = @with_kw (a = 1.5, # risk aversion - γ_x = 0.5, # production shock precision - ρ = 0.99, # correlation coefficient for θ - σ_θ = 0.5, # standard dev. of θ shock - num_firms = 100, # number of firms - σ_F = 1.5, # standard dev. of fixed costs - c = -420.0, # external opportunity cost - μ_init = 0.0, # initial value for μ - γ_init = 4.0, # initial value for γ - θ_init = 0.0, # initial value for θ - σ_x = sqrt(a / γ_x)) # standard dev. of shock +UncertaintyTrapEcon(;a = 1.5, # risk aversion + gamma_x = 0.5, # production shock precision + rho = 0.99, # correlation coefficient for theta + sigma_theta = 0.5, # standard dev. of theta shock + num_firms = 100, # number of firms + sigma_F = 1.5, # standard dev. of fixed costs + c = -420.0, # external opportunity cost + mu_init = 0.0, # initial value for mu + gamma_init = 4.0, # initial value for gamma + theta_init = 0.0, # initial value for theta + sigma_x = sqrt(a / gamma_x)) = # standard dev. of shock + (;a, gamma_x, rho, sigma_theta, num_firms, sigma_F, c, mu_init, gamma_init, theta_init, sigma_x) ``` In the results below we use this code to simulate time series for the major variables. @@ -360,17 +361,17 @@ different values of $M$ ```{code-cell} julia econ = UncertaintyTrapEcon() -@unpack ρ, σ_θ, γ_x = econ # simplify names +(;rho, sigma_theta, gamma_x) = econ # simplify names -# grid for γ and γ_{t+1} -γ = range(1e-10, 3, length = 200) +# grid for gamma and gamma_{t+1} +gamma = range(1e-10, 3, length = 200) M_range = 0:6 -γp = 1 ./ (ρ^2 ./ (γ .+ γ_x .* M_range') .+ σ_θ^2) +gammap = 1 ./ (rho^2 ./ (gamma .+ gamma_x .* M_range') .+ sigma_theta^2) labels = ["0" "1" "2" "3" "4" "5" "6"] -plot(γ, γ, lw = 2, label = "45 Degree") -plot!(γ, γp, lw = 2, label = labels) +plot(gamma, gamma, lw = 2, label = "45 Degree") +plot!(gamma, gammap, lw = 2, label = labels) plot!(xlabel = L"\gamma", ylabel = L"\gamma^\prime", legend_title = L"M", legend = :bottomright) ``` @@ -379,9 +380,9 @@ plot!(xlabel = L"\gamma", ylabel = L"\gamma^\prime", legend_title = L"M", legend tags: [remove-cell] --- @testset begin - @test γp[2,2] ≈ 0.46450522950184053 - @test γp[3,3] ≈ 0.8323524432613787 - @test γp[5,5] ≈ 1.3779664509290432 + @test gammap[2,2] ≈ 0.46450522950184053 + @test gammap[3,3] ≈ 0.8323524432613787 + @test gammap[5,5] ≈ 1.3779664509290432 end ``` @@ -396,49 +397,49 @@ is, the number of active firms and average output ```{code-cell} julia function simulate(uc, capT = 2_000) # unpack parameters - @unpack a, γ_x, ρ, σ_θ, num_firms, σ_F, c, μ_init, γ_init, θ_init, σ_x = uc + (;a, gamma_x, rho, sigma_theta, num_firms, sigma_F, c, mu_init, gamma_init, theta_init, sigma_x) = uc # draw standard normal shocks w_shocks = randn(capT) # aggregate functions - # auxiliary function ψ - function ψ(γ, μ, F) - temp1 = -a * (μ - F) - temp2 = 0.5 * a^2 / (γ + γ_x) + # auxiliary function psi + function psi(gamma, mu, F) + temp1 = -a * (mu - F) + temp2 = 0.5 * a^2 / (gamma + gamma_x) return (1 - exp(temp1 + temp2)) / a - c end # compute X, M - function gen_aggregates(γ, μ, θ) - F_vals = σ_F * randn(num_firms) - M = sum(ψ.(Ref(γ), Ref(μ), F_vals) .> 0) # counts number of active firms - if any(ψ(γ, μ, f) > 0 for f in F_vals) # ∃ an active firm - x_vals = θ .+ σ_x * randn(M) + function gen_aggregates(gamma, mu, theta) + F_vals = sigma_F * randn(num_firms) + M = sum(psi.(Ref(gamma), Ref(mu), F_vals) .> 0) # counts number of active firms + if any(psi(gamma, mu, f) > 0 for f in F_vals) # there is an active firm + x_vals = theta .+ sigma_x * randn(M) X = mean(x_vals) else X = 0.0 end - return (X = X, M = M) + return (;X, M) end # initialize dataframe - X_init, M_init = gen_aggregates(γ_init, μ_init, θ_init) - df = DataFrame(γ = γ_init, μ = μ_init, θ = θ_init, X = X_init, M = M_init) + X_init, M_init = gen_aggregates(gamma_init, mu_init, theta_init) + df = DataFrame(gamma = gamma_init, mu = mu_init, theta = theta_init, X = X_init, M = M_init) # update dataframe for t in 2:capT # unpack old variables - θ_old, γ_old, μ_old, X_old, M_old = (df.θ[end], df.γ[end], df.μ[end], df.X[end], df.M[end]) + theta_old, gamma_old, mu_old, X_old, M_old = (df.theta[end], df.gamma[end], df.mu[end], df.X[end], df.M[end]) # define new beliefs - θ = ρ * θ_old + σ_θ * w_shocks[t-1] - μ = (ρ * (γ_old * μ_old + M_old * γ_x * X_old))/(γ_old + M_old * γ_x) - γ = 1 / (ρ^2 / (γ_old + M_old * γ_x) + σ_θ^2) + theta = rho * theta_old + sigma_theta * w_shocks[t-1] + mu = (rho * (gamma_old * mu_old + M_old * gamma_x * X_old))/(gamma_old + M_old * gamma_x) + gamma = 1 / (rho^2 / (gamma_old + M_old * gamma_x) + sigma_theta^2) # compute new aggregates - X, M = gen_aggregates(γ, μ, θ) - push!(df, (γ = γ, μ = μ, θ = θ, X = X, M = M)) + X, M = gen_aggregates(gamma, mu, theta) + push!(df, (;gamma, mu, theta, X, M)) end # return @@ -459,16 +460,16 @@ Random.seed!(42); # set random seed for reproducible results ```{code-cell} julia df = simulate(econ) -plot(eachindex(df.μ), df.μ, lw = 2, label = L"\mu") -plot!(eachindex(df.θ), df.θ, lw = 2, label = L"\theta") +plot(eachindex(df.mu), df.mu, lw = 2, label = L"\mu") +plot!(eachindex(df.theta), df.theta, lw = 2, label = L"\theta") plot!(xlabel = L"x", ylabel = L"y", legend_title = "Variable", legend = :bottomright) ``` Now let's plot the whole thing together ```{code-cell} julia -len = eachindex(df.θ) -yvals = [df.θ, df.μ, df.γ, df.M] +len = eachindex(df.theta) +yvals = [df.theta, df.mu, df.gamma, df.M] vars = [L"\theta", L"\mu", L"\gamma", L"M"] plt = plot(layout = (4,1), size = (600, 600)) @@ -484,7 +485,7 @@ plot(plt) --- tags: [remove-cell] --- -mdf = DataFrame(t = eachindex(df.θ), θ = df.θ, μ = df.μ, γ = df.γ, M = df.M) +mdf = DataFrame(t = eachindex(df.theta), theta = df.theta, mu = df.mu, gamma = df.gamma, M = df.M) @testset begin @test stack(mdf, collect(2:5))[:value][3] ≈ -0.49742498224730913 atol = 1e-3 From c1f5ed2b37050d830cba31dd8a895133da73510d Mon Sep 17 00:00:00 2001 From: annabellasd Date: Tue, 10 Oct 2023 18:23:25 -0800 Subject: [PATCH 8/9] update so will merge ***ORDER matters make sure to merge continuous time last --- lectures/continuous_time/covid_sde.md | 126 ++++++++++++------------- lectures/continuous_time/seir_model.md | 82 ++++++++-------- 2 files changed, 101 insertions(+), 107 deletions(-) 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. - From 7039a567e7b9323eeab861e013ee5db0ab8f77f2 Mon Sep 17 00:00:00 2001 From: annabellasd Date: Mon, 6 Nov 2023 17:03:50 -0800 Subject: [PATCH 9/9] updates; Trang's comments --- lectures/multi_agent_models/aiyagari.md | 4 ++-- lectures/multi_agent_models/arellano.md | 12 +++++----- lectures/multi_agent_models/lake_model.md | 15 ++++++------ lectures/multi_agent_models/lucas_model.md | 2 +- lectures/multi_agent_models/markov_perf.md | 24 +++++++++---------- .../rational_expectations.md | 4 ++-- lectures/multi_agent_models/schelling.md | 1 - .../multi_agent_models/uncertainty_traps.md | 20 ++++++++-------- 8 files changed, 41 insertions(+), 41 deletions(-) diff --git a/lectures/multi_agent_models/aiyagari.md b/lectures/multi_agent_models/aiyagari.md index de0f8a5e..640ba2f2 100644 --- a/lectures/multi_agent_models/aiyagari.md +++ b/lectures/multi_agent_models/aiyagari.md @@ -226,7 +226,7 @@ Household(;r = 0.01, u = sigma == 1 ? x -> log(x) : x -> (x^(1 - sigma) - 1) / (1 - sigma), R = setup_R!(fill(-Inf, n, a_size), a_vals, s_vals, r, w, u), # -Inf is the utility of dying (0 consumption) - Q = setup_Q!(zeros(n, a_size, n), s_i_vals, z_chain)) = (; r, w, sigma, beta, z_chain, a_min, a_max, a_size, a_vals, z_size, n, s_vals, s_i_vals, u, R, Q) + Q = setup_Q!(zeros(n, a_size, n), s_i_vals, z_chain)) = (;r, w, sigma, beta, z_chain, a_min, a_max, a_size, a_vals, z_size, n, s_vals, s_i_vals, u, R, Q) function setup_Q!(Q, s_i_vals, z_chain) for next_s_i in 1:size(Q, 3) @@ -368,7 +368,7 @@ k_vals = prices_to_capital_stock.(Ref(am), r_vals) # Plot against demand for capital by firms demand = rd.(k_vals) -labels = ["demand for capital" "supply of capital"] +labels = ["demand for capital" "supply of capital"] plot(k_vals, [demand r_vals], label = labels, lw = 2, alpha = 0.6) plot!(xlabel = "capital", ylabel = "interest rate", xlim = (2, 14), ylim = (0.0, 0.1)) ``` diff --git a/lectures/multi_agent_models/arellano.md b/lectures/multi_agent_models/arellano.md index e9e92a40..6cc33a9f 100644 --- a/lectures/multi_agent_models/arellano.md +++ b/lectures/multi_agent_models/arellano.md @@ -590,14 +590,14 @@ using DataFrames, Plots Compute the value function, policy and equilibrium prices ```{code-cell} julia -ae = ArellanoEconomy(beta = .953, # time discount rate +ae = ArellanoEconomy(beta = .953, # time discount rate gamma = 2., # risk aversion - r = 0.017, # international interest rate - rho = .945, # persistence in output - eta = 0.025, # st dev of output shock + r = 0.017, # international interest rate + rho = .945, # persistence in output + eta = 0.025, # st dev of output shock theta = 0.282, # prob of regaining access - ny = 21, # number of points in y grid - nB = 251) # number of points in B grid + ny = 21, # number of points in y grid + nB = 251) # number of points in B grid # now solve the model on the grid. vfi!(ae) diff --git a/lectures/multi_agent_models/lake_model.md b/lectures/multi_agent_models/lake_model.md index 1866c63e..63e46e5b 100644 --- a/lectures/multi_agent_models/lake_model.md +++ b/lectures/multi_agent_models/lake_model.md @@ -221,7 +221,7 @@ function transition_matrices(lm) A = [(1 - lambda) * (1 - d) + b (1 - d) * alpha + b (1 - d) * lambda (1 - d) * (1 - alpha)] A_hat = A ./ (1 + g) - return (A, A_hat) + return (;A, A_hat) end function rate_steady_state(lm) @@ -630,7 +630,7 @@ Following {cite}`davis2006flow`, we set $\alpha$, the hazard rate of leaving emp We will make use of (with some tweaks) the code we wrote in the {doc}`McCall model lecture <../dynamic_programming/mccall_model>`, embedded below for convenience. ```{code-cell} julia -function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1e-5,iter = 2_000) +function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1e-5, iter = 2_000) (;alpha, beta, sigma, c, gamma, w, E, u) = mcm # necessary objects @@ -671,13 +671,14 @@ u(c, sigma) = c > 0 ? (c^(1 - sigma) - 1) / (1 - sigma) : -10e-6 # model constructor McCallModel(;alpha = 0.2, - beta = 0.98, # discount rate + beta = 0.98, # discount rate gamma = 0.7, - c = 6.0, # unemployment compensation + c = 6.0, # unemployment compensation sigma = 2.0, - u = u, # utility function - w = range(10, 20, length = 60), # w= wage values; E= distribution over wage values - E = Expectation(BetaBinomial(59, 600, 400))) = (;alpha, beta, gamma, c, sigma, u, w, E) + u = u, # utility function + w = range(10, 20, length = 60), # w= wage values + E = Expectation(BetaBinomial(59, 600, 400))) = # E = distribution over wage values + (;alpha, beta, gamma, c, sigma, u, w, E) ``` diff --git a/lectures/multi_agent_models/lucas_model.md b/lectures/multi_agent_models/lucas_model.md index 2b26c385..c3031492 100644 --- a/lectures/multi_agent_models/lucas_model.md +++ b/lectures/multi_agent_models/lucas_model.md @@ -413,7 +413,7 @@ function LucasTree(;gamma = 2.0, h[i] = beta * mean((y^alpha .* shocks).^(1 - gamma)) end - return (gamma = gamma, beta = beta, alpha = alpha, sigma = sigma, phi = phi, grid = grid, shocks = shocks, h = h) + return (;gamma, beta, alpha, sigma, phi, grid, shocks, h) end # approximate Lucas operator, which returns the updated function Tf on the grid diff --git a/lectures/multi_agent_models/markov_perf.md b/lectures/multi_agent_models/markov_perf.md index 06fcfcf9..24fcb846 100644 --- a/lectures/multi_agent_models/markov_perf.md +++ b/lectures/multi_agent_models/markov_perf.md @@ -440,13 +440,13 @@ A = I + zeros(3, 3) B1 = [0.0, 1.0, 0.0] B2 = [0.0, 0.0, 1.0] -R1 = [ 0.0 -a0 / 2.0 0.0; - -a0 / 2.0 a1 a1 / 2.0; - 0.0 a1 / 2.0 0.0] +R1 = [ 0.0 -a0 / 2.0 0.0; + -a0 / 2.0 a1 a1 / 2.0; + 0.0 a1 / 2.0 0.0 ] -R2 = [ 0.0 0.0 -a0 / 2.0; - 0.0 0.0 a1 / 2.0; - -a0 / 2.0 a1 / 2.0 a1] +R2 = [ 0.0 0.0 -a0 / 2.0; + 0.0 0.0 a1 / 2.0; + -a0 / 2.0 a1 / 2.0 a1 ] Q1 = Q2 = gamma S1 = S2 = W1 = W2 = M1 = M2 = 0.0 @@ -480,8 +480,8 @@ In particular, let's take F2 as computed above, plug it into {eq}`eq_mpe_p1p` an We hope that the resulting policy will agree with F1 as computed above ```{code-cell} julia -Lamda1 = A - (B2 * F2) -lq1 = QuantEcon.LQ(Q1, R1, Lamda1, B1, bet=beta) +Lambda1 = A - (B2 * F2) +lq1 = QuantEcon.LQ(Q1, R1, Lambda1, B1, bet=beta) P1_ih, F1_ih, d = stationary_values(lq1) F1_ih ``` @@ -493,7 +493,7 @@ tags: [remove-cell] @testset begin @test P1_ih[2, 2] ≈ 5.441368459897164 @test d ≈ 0.0 - @test Lamda1[1, 1] ≈ 1.0 && Lamda1[3, 2] ≈ -0.07584666305807419 + @test Lambda1[1, 1] ≈ 1.0 && Lambda1[3, 2] ≈ -0.07584666305807419 @test F1_ih ≈ [-0.6684661291052371 0.29512481789806305 0.07584666292394007] @test isapprox(F1, F1_ih, atol=1e-7) # Make sure the test below comes up true. end @@ -860,9 +860,9 @@ we set up the matrices as follows: ```{code-cell} julia # create matrices needed to compute the Nash feedback equilibrium -A = [delta_1 0 -delta_1 * b[1]; - 0 delta_1 -delta_1 * b[2]; - 0 0 1] +A = [delta_1 0 -delta_1 * b[1]; + 0 delta_1 -delta_1 * b[2]; + 0 0 1 ] B1 = delta_1 * [1 -D[1, 1]; 0 -D[2, 1]; diff --git a/lectures/multi_agent_models/rational_expectations.md b/lectures/multi_agent_models/rational_expectations.md index 2215b90e..a927c0ef 100644 --- a/lectures/multi_agent_models/rational_expectations.md +++ b/lectures/multi_agent_models/rational_expectations.md @@ -670,9 +670,9 @@ kappa0 = 95.5 kappa1 = 0.95 # formulate the LQ problem -A = [1 0 0 +A = [1 0 0 0 kappa1 kappa0 - 0 0 1] + 0 0 1 ] B = [1.0, 0.0, 0.0] diff --git a/lectures/multi_agent_models/schelling.md b/lectures/multi_agent_models/schelling.md index 4bb0f874..63bbcd37 100644 --- a/lectures/multi_agent_models/schelling.md +++ b/lectures/multi_agent_models/schelling.md @@ -8,7 +8,6 @@ kernelspec: language: julia name: julia-1.9 --- -##AT THIS ONE BUT LUCAS (schelling)= ```{raw} html
diff --git a/lectures/multi_agent_models/uncertainty_traps.md b/lectures/multi_agent_models/uncertainty_traps.md index 477ef665..a0a531a1 100644 --- a/lectures/multi_agent_models/uncertainty_traps.md +++ b/lectures/multi_agent_models/uncertainty_traps.md @@ -236,16 +236,16 @@ using DataFrames, LaTeXStrings, Parameters, Plots ``` ```{code-cell} julia -UncertaintyTrapEcon(;a = 1.5, # risk aversion - gamma_x = 0.5, # production shock precision - rho = 0.99, # correlation coefficient for theta - sigma_theta = 0.5, # standard dev. of theta shock - num_firms = 100, # number of firms - sigma_F = 1.5, # standard dev. of fixed costs - c = -420.0, # external opportunity cost - mu_init = 0.0, # initial value for mu - gamma_init = 4.0, # initial value for gamma - theta_init = 0.0, # initial value for theta +UncertaintyTrapEcon(;a = 1.5, # risk aversion + gamma_x = 0.5, # production shock precision + rho = 0.99, # correlation coefficient for theta + sigma_theta = 0.5, # standard dev. of theta shock + num_firms = 100, # number of firms + sigma_F = 1.5, # standard dev. of fixed costs + c = -420.0, # external opportunity cost + mu_init = 0.0, # initial value for mu + gamma_init = 4.0, # initial value for gamma + theta_init = 0.0, # initial value for theta sigma_x = sqrt(a / gamma_x)) = # standard dev. of shock (;a, gamma_x, rho, sigma_theta, num_firms, sigma_F, c, mu_init, gamma_init, theta_init, sigma_x) ```