134134
135135Interpretation:
136136
137- * draw $q$ from a β distribution with shape parameters $(a, b)$
137+ * draw $q$ from a $\beta$ distribution with shape parameters $(a, b)$
138138* run $n$ independent binary trials, each with success probability $q$
139139* $p(k \, |\, n, a, b)$ is the probability of $k$ successes in these $n$ trials
140140
@@ -156,7 +156,6 @@ using Test
156156``` {code-cell} julia
157157using LaTeXStrings, Plots, QuantEcon, Distributions
158158
159-
160159n = 50
161160a_vals = [0.5, 1, 100]
162161b_vals = [0.5, 1, 100]
@@ -182,41 +181,37 @@ Implementation:
182181The code for solving the DP problem described above is found below:
183182
184183``` {code-cell} julia
185- function CareerWorkerProblem(;β = 0.95,
184+ function CareerWorkerProblem(; beta = 0.95,
186185 B = 5.0,
187186 N = 50,
188187 F_a = 1.0,
189188 F_b = 1.0,
190189 G_a = 1.0,
191190 G_b = 1.0)
192- θ = range(0, B, length = N)
193- ϵ = copy(θ )
194- dist_F = BetaBinomial(N- 1, F_a, F_b)
195- dist_G = BetaBinomial(N- 1, G_a, G_b)
191+ theta = range(0, B, length = N)
192+ epsilon = copy(theta )
193+ dist_F = BetaBinomial(N - 1, F_a, F_b)
194+ dist_G = BetaBinomial(N - 1, G_a, G_b)
196195 F_probs = pdf.(dist_F, support(dist_F))
197196 G_probs = pdf.(dist_G, support(dist_G))
198- F_mean = sum(θ .* F_probs)
199- G_mean = sum(ϵ .* G_probs)
200- return (β = β, N = N, B = B, θ = θ, ϵ = ϵ,
201- F_probs = F_probs, G_probs = G_probs,
202- F_mean = F_mean, G_mean = G_mean)
197+ F_mean = sum(theta .* F_probs)
198+ G_mean = sum(epsilon .* G_probs)
199+ return (; beta, N, B, theta, epsilon, F_probs, G_probs, F_mean, G_mean)
203200end
204201
205202function update_bellman!(cp, v, out; ret_policy = false)
206203
207204 # new life. This is a function of the distribution parameters and is
208205 # always constant. No need to recompute it in the loop
209- v3 = (cp.G_mean + cp.F_mean .+ cp.β .*
210- cp.F_probs' * v * cp.G_probs)[1] # do not need 1 element array
206+ v3 = (cp.G_mean + cp.F_mean .+ cp.beta .* cp.F_probs' * v * cp.G_probs)[1] # do not need 1 element array
211207
212- for j in 1:cp.N
213- for i in 1:cp.N
208+ for j in 1:( cp.N)
209+ for i in 1:( cp.N)
214210 # stay put
215- v1 = cp.θ [i] + cp.ϵ [j] + cp.β * v[i, j]
211+ v1 = cp.theta [i] + cp.epsilon [j] + cp.beta * v[i, j]
216212
217213 # new job
218- v2 = (cp.θ[i] .+ cp.G_mean .+ cp.β .*
219- v[i, :]' * cp.G_probs)[1] # do not need a single element array
214+ v2 = (cp.theta[i] .+ cp.G_mean .+ cp.beta .* v[i, :]' * cp.G_probs)[1] # do not need a single element array
220215
221216 if ret_policy
222217 if v1 > max(v2, v3)
@@ -234,7 +229,6 @@ function update_bellman!(cp, v, out; ret_policy = false)
234229 end
235230end
236231
237-
238232function update_bellman(cp, v; ret_policy = false)
239233 out = similar(v)
240234 update_bellman!(cp, v, out, ret_policy = ret_policy)
@@ -274,8 +268,9 @@ v_init = fill(100.0, wp.N, wp.N)
274268func(x) = update_bellman(wp, x)
275269v = compute_fixed_point(func, v_init, max_iter = 500, verbose = false)
276270
277- plot(linetype = :surface, wp.θ, wp.ϵ, transpose(v), xlabel=L"\theta", ylabel=L"\epsilon",
278- seriescolor=:plasma, gridalpha = 1)
271+ plot(linetype = :surface, wp.theta, wp.epsilon, transpose(v), xlabel = L"\theta",
272+ ylabel = L"\epsilon",
273+ seriescolor = :plasma, gridalpha = 1)
279274```
280275
281276The optimal policy can be represented as follows (see {ref}` Exercise 3 <career_ex3> ` for code).
@@ -373,31 +368,31 @@ G = DiscreteRV(wp.G_probs)
373368
374369function gen_path(T = 20)
375370 i = j = 1
376- θ_ind = Int[]
377- ϵ_ind = Int[]
371+ theta_ind = Int[]
372+ epsilon_ind = Int[]
378373
379- for t= 1:T
374+ for t in 1:T
380375 # do nothing if stay put
381376 if optimal_policy[i, j] == 2 # new job
382377 j = rand(G)[1]
383378 elseif optimal_policy[i, j] == 3 # new life
384379 i, j = rand(F)[1], rand(G)[1]
385380 end
386- push!(θ_ind , i)
387- push!(ϵ_ind , j)
381+ push!(theta_ind , i)
382+ push!(epsilon_ind , j)
388383 end
389- return wp.θ[θ_ind ], wp.ϵ[ϵ_ind ]
384+ return wp.theta[theta_ind ], wp.epsilon[epsilon_ind ]
390385end
391386
392387plot_array = Any[]
393388for i in 1:2
394- θ_path, ϵ_path = gen_path()
395- plt = plot(ϵ_path , label= L"\epsilon")
396- plot!(plt, θ_path , label= L"\theta")
397- plot!(plt, legend= :bottomright)
389+ theta_path, epsilon_path = gen_path()
390+ plt = plot(epsilon_path , label = L"\epsilon")
391+ plot!(plt, theta_path , label = L"\theta")
392+ plot!(plt, legend = :bottomright)
398393 push!(plot_array, plt)
399394end
400- plot(plot_array..., layout = (2,1))
395+ plot(plot_array..., layout = (2, 1))
401396```
402397
403398``` {code-cell} julia
@@ -430,7 +425,6 @@ function gen_first_passage_time(optimal_policy)
430425 end
431426end
432427
433-
434428M = 25000
435429samples = zeros(M)
436430for i in 1:M
@@ -448,12 +442,12 @@ tags: [remove-cell]
448442end
449443```
450444
451- To compute the median with $\beta=0.99$ instead of the default value $\beta=0.95$, replace ` wp=CareerWorkerProblem() ` with ` wp=CareerWorkerProblem(β =0.99) ` .
445+ To compute the median with $\beta=0.99$ instead of the default value $\beta=0.95$, replace ` wp=CareerWorkerProblem() ` with ` wp=CareerWorkerProblem(beta =0.99) ` .
452446
453447The medians are subject to randomness, but should be about 7 and 14 respectively. Not surprisingly, more patient workers will wait longer to settle down to their final job.
454448
455449``` {code-cell} julia
456- wp2 = CareerWorkerProblem(β= 0.99)
450+ wp2 = CareerWorkerProblem(beta = 0.99)
457451
458452v2, optimal_policy2 = solve_wp(wp2)
459453
@@ -485,30 +479,30 @@ lvls = [0.5, 1.5, 2.5, 3.5]
485479x_grid = range(0, 5, length = 50)
486480y_grid = range(0, 5, length = 50)
487481
488- contour(x_grid, y_grid, optimal_policy', fill= true, levels= lvls,color = :Blues,
489- fillalpha= 1, cbar = false)
490- contour!(xlabel= L"\theta", ylabel= L"\epsilon")
491- annotate!([(1.8,2.5, text("new life", 14, :white, :center))])
492- annotate!([(4.5,2.5, text("new job", 14, :center))])
493- annotate!([(4.0,4.5, text("stay put", 14, :center))])
482+ contour(x_grid, y_grid, optimal_policy', fill = true, levels = lvls, color = :Blues,
483+ fillalpha = 1, cbar = false)
484+ contour!(xlabel = L"\theta", ylabel = L"\epsilon")
485+ annotate!([(1.8, 2.5, text("new life", 14, :white, :center))])
486+ annotate!([(4.5, 2.5, text("new job", 14, :center))])
487+ annotate!([(4.0, 4.5, text("stay put", 14, :center))])
494488```
495489
496490Now, we need only swap out for the new parameters
497491
498492``` {code-cell} julia
499- wp = CareerWorkerProblem(G_a= 100.0, G_b= 100.0); # use new params
493+ wp = CareerWorkerProblem(G_a = 100.0, G_b = 100.0); # use new params
500494v, optimal_policy = solve_wp(wp)
501495
502496lvls = [0.5, 1.5, 2.5, 3.5]
503497x_grid = range(0, 5, length = 50)
504498y_grid = range(0, 5, length = 50)
505499
506- contour(x_grid, y_grid, optimal_policy', fill= true, levels= lvls,color = :Blues,
507- fillalpha= 1, cbar = false)
508- contour!(xlabel= L"\theta", ylabel= L"\epsilon")
509- annotate!([(1.8,2.5, text("new life", 14, :white, :center))])
510- annotate!([(4.5,2.5, text("new job", 14, :center))])
511- annotate!([(4.0,4.5, text("stay put", 14, :center))])
500+ contour(x_grid, y_grid, optimal_policy', fill = true, levels = lvls, color = :Blues,
501+ fillalpha = 1, cbar = false)
502+ contour!(xlabel = L"\theta", ylabel = L"\epsilon")
503+ annotate!([(1.8, 2.5, text("new life", 14, :white, :center))])
504+ annotate!([(4.5, 2.5, text("new job", 14, :center))])
505+ annotate!([(4.0, 4.5, text("stay put", 14, :center))])
512506```
513507
514508You will see that the region for which the worker
0 commit comments