using QuantEcon import QuantEcon: solve using Distributions using PyPlot using Roots type JobSearchModel # Parameters w::Vector{Float64} w_pdf::Vector{Float64} beta::Float64 alpha::Float64 gamma::Float64 rho::Float64 # Internal variables u::Function ddp::DiscreteDP num_states::Integer num_actions::Integer rej::Integer acc::Integer function JobSearchModel(w::Vector, w_pdf::Vector, beta::Float64; alpha::Float64=0., gamma::Float64=1., rho::Float64=0.) # Utility function function u(y::Vector{Float64}) small_number = -9999999 nonpositive = (y .<= 0) if rho == 1 util = log(y) else util = (y.^(1 - rho) - 1)/(1 - rho) end util[nonpositive] = small_number return util end u_w = u(w) num_states = length(w) + 1 num_actions = 2 rej, acc = 1, 2 # Reward array R0 = zeros(num_states, num_actions) R0[1:end-1, acc] = u_w # Transition probability array Q = zeros(num_states, num_actions, num_states) # Reject for s in 1:num_states Q[s, rej, 1:end-1] = w_pdf * gamma Q[s, rej, end] = 1 - gamma end # Accept for s in 1:num_states-1 Q[s, acc, s] = 1 - alpha Q[s, acc, end] = alpha end Q[end, acc, 1:end-1] = w_pdf ddp = DiscreteDP(R0, Q, beta) js = new(w, w_pdf, beta, alpha, gamma, rho, u, ddp, num_states, num_actions, rej, acc) return js end end; function solve(js::JobSearchModel, c::Float64) n, m = js.num_states, js.num_actions rej = js.rej js.ddp.R[:, rej] = js.u([c])[1] js.ddp.R[end, :] = js.u([c])[1] res = solve(js.ddp, PFI) V::Vector{Float64} = res.v[1:end-1] # Values of jobs U::Float64 = res.v[end] # Value of unemployed C::Vector{Int64} = res.sigma[1:end-1] - 1 # Optimal policy lamb::Float64 = dot(js.w_pdf, C) * js.gamma pi::Vector{Float64} = [js.alpha, lamb] pi /= sum(pi) # Stationary distribution return V, U, C, pi end; w = linspace(0, 175, 201) # wage grid # compute probability of each wage level logw_dist = Normal(log(20), 1) logw_dist_cdf = cdf(logw_dist, log(w)) logw_dist_pdf = logw_dist_cdf[2:end] - logw_dist_cdf[1:end-1] logw_dist_pdf /= sum(logw_dist_pdf) w = Array((w[2:end] + w[1:end-1])/2); gamma = 1. alpha = 0.013 # Monthly alpha_q = (1-(1-alpha)^3) # Quarterly beta = 0.99 rho = 2. # risk-aversion; js = JobSearchModel(w, logw_dist_pdf, beta, alpha=alpha_q, gamma=gamma, rho=rho); c = 40. V, U, C, pi = solve(js, c); s_star = length(w) - sum(C) + 1 println("Optimal policy: Accept if and only if w >= $(w[s_star])") fig, ax = subplots(figsize=(6, 4)) ax[:plot](w, V, label=L"$V$") ax[:plot]((w[1], w[end]), (U, U), "r--", label=L"$U$") ax[:set_title]("Optimal value function") ax[:set_xlabel]("Wage") ax[:set_ylabel]("Value") legend(loc=2) show() type UnemploymentInsurancePolicy w::Vector{Float64} w_pdf::Vector{Float64} beta::Float64 alpha::Float64 gamma::Float64 rho::Float64 function UnemploymentInsurancePolicy(w::Vector, w_pdf::Vector, beta::Float64; alpha::Float64=0., gamma::Float64=1., rho::Float64=0.) uip = new(w, w_pdf, beta, alpha_q, gamma, rho) return uip end end; function solve_job_search_model(uip::UnemploymentInsurancePolicy, c::Float64, T::Float64) js = JobSearchModel(uip.w-T, uip.w_pdf, uip.beta, alpha=uip.alpha, gamma=uip.gamma, rho=uip.rho) V, U, C, pi = solve(js, c-T) return V, U, C, pi end; function budget_balance(uip::UnemploymentInsurancePolicy, c::Float64, T::Float64) V, U, C, pi = solve_job_search_model(uip, c, T) return T - pi[1]*c end; function implement(uip::UnemploymentInsurancePolicy, c::Float64) # Budget balancing tax given c T = fzero(T -> budget_balance(uip, c, T), 0., c, xtolrel=1e-3) V, U, C, pi = solve_job_search_model(uip, c, T) EV = dot(C .* V, uip.w_pdf) / (dot(C, uip.w_pdf)) W = pi[1] * U + pi[2] * EV return T, W, pi end; uip = UnemploymentInsurancePolicy(w, logw_dist_pdf, beta, alpha=alpha_q, gamma=gamma, rho=rho); grid_size = 26 cvec = linspace(5, 135, grid_size) Ts, Ws = Array(Float64, grid_size), Array(Float64, grid_size) pis = Array(Float64, 2, grid_size) for (i, c) in enumerate(cvec) T, W, pi = implement(uip, c) Ts[i], Ws[i], pis[:, i] = T, W, pi end i_max = indmax(Ws) println("Optimal unemployment benefit: $(cvec[i_max])") function plot(ax, y_vec, title) ax[:plot](cvec, y_vec) ax[:set_xlabel](L"$c$") ax[:vlines](cvec[i_max], ax[:get_ylim]()[1], y_vec[i_max], "k", "-.") ax[:set_title](title) end fig, axes = subplots(2, 2) plot(axes[1, 1], Ws, "Welfare") plot(axes[1, 2], Ts, "Taxes") plot(axes[2, 1], vec(pis[2, :]), "Employment Rate") plot(axes[2, 2], vec(pis[1, :]), "Unemployment Rate") tight_layout() show() fig, ax = subplots(figsize=(6, 4)) plot(ax, cvec-Ts, "Net Compensation") show()