Cost of posting a vacancy: $C_v$; prob. of filling a vacancy: $q$; cost per worker $C_v/q$
workers' outside option: $b$ unemployment benefit
workers' utility cost of working $h$ hours: $\xi_0 (h)^{\xi_1}$, $\xi_1>0, \xi_1>1$
labour market: $n=M(U,V)$, $n$=matches, $U$ job searchers, $V$ vacancies
In [119]:
    
using NLsolve
using PyPlot
    
In [120]:
    
# model parameters
Cv=.1;
psi0=.025;
psi1=1.05;
alpha=0.4;
bbeta=0.2;
epsilon=1.0;
gamma=.7;
mu=0.7;
z=1.0;
b=.1;
g=0.1;
    
In [121]:
    
function q(mm::Float64,gg::Float64,n)
    (mm*((1-n)./n).^gg).^(1/(1-gg))
end
    
    Out[121]:
In [122]:
    
function ztilda(zed::Float64,zedlow::Float64,urate)
    ztil=(1-urate).*zed+urate.*zedlow
end
    
    Out[122]:
In [123]:
    
function ipol(nbar::Float64,n::Float64)
    if nbar-n>0
        i=0;
        else 
        i=0;
    end
end
    
    Out[123]:
In [124]:
    
function twosectorsbaseline!(x,residuals)
    # x[1]=n1; x[2]=h1; x[3]=n2; x[4]=h2; x[5]=q; 
    
    # sector 1:
    # n1:
    #residuals[1]=alpha*epsilon*z*x[2]*(z*x[1]*x[2])^(alpha-1)-(b+psi0*x[2]^(1+psi1))-Cv/x[5]
    residuals[1]=epsilon*alpha*bbeta*(1/x[1])*((x[1]^bbeta)*x[2]^(1-bbeta))^(alpha)-(b+psi0*x[2]^(1+psi1))-Cv/x[5]
    # h1:
    #residuals[2]=alpha*epsilon*z*(z*x[1]*x[2])^(alpha-1)-(b/x[2]+psi0*x[2]^(psi1))-(-b/x[2]+psi0*psi1*x[2]^(psi1))
    residuals[2]=epsilon*alpha*(1-bbeta)*(1/x[2])*((x[1]^bbeta)*x[2]^(1-bbeta))^(alpha)-(b/x[2]+psi0*x[2]^(psi1))-(-b/x[2]+psi0*psi1*x[2]^(psi1))
    # sector 2:    
    # n2:
    residuals[3]=1*alpha*bbeta*(1/x[3])*((x[3]^bbeta)*x[4]^(1-bbeta))^(alpha)-(b+psi0*x[4]^(1+psi1))-Cv/x[5]    
    # h2:
    residuals[4]=1*alpha*(1-bbeta)*(1/x[4])*((x[3]^bbeta)*x[4]^(1-bbeta))^(alpha)-(b/x[4]+psi0*x[4]^(psi1))-(-b/x[4]+psi0*psi1*x[4]^(psi1))
    
    # q:
    residuals[5]=x[5]-(mu*((1-x[1]-x[3])./(x[1]+x[3])).^gamma).^(1/(1-gamma))
    
end
    
    Out[124]:
In [125]:
    
n1baseline=Array{Float64,1}(0)
n2baseline=Array{Float64,1}(0)
h1baseline=Array{Float64,1}(0)
h2baseline=Array{Float64,1}(0)
qbaseline=Array{Float64,1}(0)
for eps=1.0:-.01:0.8
    epsilon=eps;
    r=nlsolve(twosectorsbaseline!, [ .2; .2; .5; .5; .1])
    push!(n1baseline,r.zero[1])
    push!(h1baseline,r.zero[2])
    push!(n2baseline,r.zero[3])
    push!(h2baseline,r.zero[4])
    push!(qbaseline,r.zero[5])
end
    
In [126]:
    
hbar=h1baseline[1]
nbar=n1baseline[1]
function twosectorspolicy!(x,residuals)
     # x[1]=n1; x[2]=h1; x[3]=n2; x[4]=h2; x[5]=q; 
    
    # sector 1:
    # n1:
    residuals[1]=epsilon*alpha*bbeta*(1/x[1])*((x[1]^bbeta)*x[2]^(1-bbeta))^(alpha)-(b-g*(hbar-x[2])+psi0*x[2]^(1+psi1))-Cv/x[5]
    # h1:
    residuals[2]=epsilon*alpha*(1-bbeta)*(1/x[2])*((x[1]^bbeta)*x[2]^(1-bbeta))^(alpha)-(b-g*(hbar-x[2])+psi0*x[2]^(1+psi1))/x[2]-(-b/x[2]+g*hbar/x[2]+psi0*psi1*x[2]^(psi1))
    # sector 2:    
    # n2:
    residuals[3]=1*alpha*bbeta*(1/x[3])*((x[3]^bbeta)*x[4]^(1-bbeta))^(alpha)-(b-g*(hbar-x[4])+psi0*x[4]^(1+psi1))-Cv/x[5]    
    # h2:
    residuals[4]=1*alpha*(1-bbeta)*(1/x[4])*((x[3]^bbeta)*x[4]^(1-bbeta))^(alpha)-(b-g*(hbar-x[4])+psi0*x[4]^(1+psi1))/x[4]-(-b/x[4]+g*hbar/x[4]+psi0*psi1*x[4]^(psi1))
    
    # q:
    residuals[5]=x[5]-(mu*((1-x[1]-x[3])./(x[1]+x[3])).^gamma).^(1/(1-gamma))
    
end
    
    Out[126]:
In [127]:
    
n1policy=Array{Float64,1}(0)
n2policy=Array{Float64,1}(0)
h1policy=Array{Float64,1}(0)
h2policy=Array{Float64,1}(0)
qpolicy=Array{Float64,1}(0)
for eps=1.0:-.01:0.8
    epsilon=eps;
    #r=nlsolve(twosectorspolicy!, [ n1baseline[end] ; h1baseline[end] ; n2baseline[end]; h2baseline[end]; qbaseline[end]])
    r=nlsolve(twosectorspolicy!, [ .2; .2; .5; .5; .1])
    push!(n1policy,r.zero[1])
    push!(h1policy,r.zero[2])
    push!(n2policy,r.zero[3])
    push!(h2policy,r.zero[4])
    push!(qpolicy,r.zero[5])
end
    
In [ ]:
    
    
In [ ]:
    
    
In [ ]:
    
    
In [ ]:
    
    
In [128]:
    
figure("employmentsteadystate",figsize=(8,6))
fig, ax = subplots()
ax[:plot](1,marker="o", linewidth=1.0, color="red", alpha=0.9, label="Sector 1")
ax[:plot](1, marker="o", linewidth=1.0, color="blue", alpha=0.5, label="Sector 2")
ax[:legend](loc="lower left")
ax[:set_title]("Employment by sector")
ax[:set_xticks]([0, .1, .2])
#ax[:set_yticks]([0.7,0.8,0.9,1.0, 1.2])
xlabel(L"$\Delta_\epsilon$")
savefig("employmentsteadystate.png")
close("employmentsteadystate")
    
    
In [129]:
    
figure("employmentbaseline",figsize=(8,6))
delta=1-collect(1.0:-.01:0.8)
fig, ax = subplots()
ax[:plot](delta,n1baseline./n1baseline[1], linewidth=2, color="red", alpha=0.9, label="Sector 1")
ax[:plot](delta,n2baseline./n2baseline[1], linewidth=2, color="blue", alpha=0.9, label="Sector 2")
ax[:legend](loc="lower left")
ax[:set_title]("Employment by sector")
ax[:set_xticks]([0, .1, .2])
#ax[:set_yticks]([0.7,0.8,0.9,1.0, 1.2])
xlabel(L"$\Delta_\epsilon$")
savefig("employmentbaseline.png")
close("employmentbaseline")
    
    
In [130]:
    
figure("hoursbaseline",figsize=(8,6))
fig, ax = subplots()
ax[:plot](delta,h1baseline./h1baseline[1], linewidth=2, color="red", alpha=0.9, label="Sector 1")
ax[:plot](delta,h2baseline./h2baseline[1], linewidth=2, color="blue", alpha=0.9, label="Sector 2")
ax[:legend](loc="lower left")
ax[:set_title]("Hours")
ax[:set_xticks]([0, .1, .2])
xlabel(L"$\Delta_\epsilon$")
savefig("hoursbaseline.png")
close("hoursbaseline")
    
    
In [131]:
    
figure("employmentpolicy",figsize=(8,6))
fig, ax = subplots()
ax[:plot](delta,n1baseline./n1baseline[1], linewidth=2, color="red", alpha=0.9, label="Sector 1")
ax[:plot](delta,n2baseline./n2baseline[1], linewidth=2, color="blue", alpha=0.9, label="Sector 2")
ax[:plot](delta,n1policy./n1policy[1], linewidth=2, linestyle="--", color="red", alpha=0.5, label="Sector 1, policy")
ax[:plot](delta,n2policy./n2policy[1], linewidth=2, linestyle="--", color="blue", alpha=0.5, label="Sector 2, policy")
ax[:legend](loc="lower left")
ax[:set_title]("Employment")
ax[:set_xticks]([0, .1, .2])
xlabel(L"$\Delta_\epsilon$")
savefig("employmentpolicy.png")
close("employmentpolicy")
    
    
In [132]:
    
figure("hourspolicy",figsize=(8,6))
fig, ax = subplots()
ax[:plot](delta,h1baseline./h1baseline[1], linewidth=2, color="red", alpha=0.9, label="Sector 1")
ax[:plot](delta,h2baseline./h2baseline[1], linewidth=2, color="blue", alpha=0.9, label="Sector 2")
ax[:plot](delta,h1policy./h1policy[1], linewidth=2, linestyle="--", color="red", alpha=0.5, label="Sector 1, policy")
ax[:plot](delta,h2policy./h2policy[1], linewidth=2, linestyle="--", color="blue", alpha=0.5, label="Sector 2, policy")
ax[:legend](loc="lower left")
ax[:set_title]("Hours")
ax[:set_xticks]([0, .1, .2])
xlabel(L"$\Delta_\epsilon$")
savefig("hourspolicy.png")
close("hourspolicy")
    
    
In [133]:
    
figure("aggregateemployment",figsize=(8,6))
fig, ax = subplots()
ax[:plot](delta,(n1baseline+n2baseline)./(n1baseline[1]+n2baseline[1]), linewidth=2, color="black", alpha=0.9, label="baseline")
ax[:plot](delta,(n1policy+n2policy)./(n1policy[1]+n2policy[1]), linewidth=2, linestyle="--", color="black", alpha=0.5, label="policy")
ax[:legend](loc="upper right")
ax[:set_title]("Aggregate employment")
ax[:set_xticks]([0, .1, .2])
xlabel(L"$\Delta_\epsilon$")
savefig("aggregateemployment.png")
close("aggregateemployment")
    
    
In [134]:
    
figure("aggregatehours",figsize=(8,6))
fig, ax = subplots()
ax[:plot](delta,(h1baseline+n2baseline)./(h1baseline[1]+n2baseline[1]), linewidth=2, color="black", alpha=0.9, label="baseline")
ax[:plot](delta,(h1policy+n2policy)./(h1policy[1]+n2policy[1]), linewidth=2, linestyle="--", color="black", alpha=0.5, label="policy")
ax[:legend](loc="upper right")
ax[:set_title]("Aggregate hours")
ax[:set_xticks]([0, .1, .2])
xlabel(L"$\Delta_\epsilon$")
savefig("aggregatehours.png")
close("aggregatehours")
    
    
In [135]:
    
outputbaseline=collect(1.0:-.01:0.8).*((n1baseline.^bbeta).*h1baseline.^(1-bbeta)).^alpha+
((n2baseline.^bbeta).*h2baseline.^(1-bbeta)).^alpha;
outputpolicy=collect(1.0:-.01:0.8).*((n1policy.^bbeta).*h1policy.^(1-bbeta)).^alpha+
((n2policy.^bbeta).*h2policy.^(1-bbeta)).^alpha;
    
In [136]:
    
figure("output",figsize=(8,6))
fig, ax = subplots()
ax[:plot](delta,outputbaseline./outputbaseline[1], linewidth=2, color="black", alpha=0.9, label="baseline")
ax[:plot](delta,outputpolicy./outputpolicy[1], linewidth=2, linestyle="--", color="black", alpha=0.5, label="policy")
ax[:legend](loc="upper right")
ax[:set_title]("Output")
ax[:set_xticks]([0, .1, .2])
xlabel(L"$\Delta_\epsilon$")
savefig("output.png")
close("output")