In [131]:
using NLsolve
using PyPlot

In [132]:
# model parameters
Cv=0.2;
psi1=1.01;
alpha=0.5;
epsilon=1.0;
hbar=1.0;
gamma=.85;
mu=0.9;
z=1.0;
zlow=.5;

In [133]:
function q(mm::Float64,gg::Float64,n)
    (mm*((1-n)./n).^gg).^(1/(1-gg))
end


Out[133]:
q (generic function with 1 method)

In [145]:
function ztilda(zed::Float64,zedlow::Float64,urate)
    ztil=(1-urate).*zed+urate.*zedlow
end


Out[145]:
ztilda (generic function with 2 methods)

In [135]:
function twoperiodeqflexible!(x,residuals)
    # x[1]=n1; x[2]=h1; x[3]=n2z; x[4]=h2z; x[5]=n2zlow; x[6]=h2zlow; x[7]=q1; x[8]=q2z; 
    
    # Last period FOC:
    
    # n2:
    residuals[1]=alpha*1*ztilda(z,zlow,1-x[1])*x[4]*(z*x[3]*x[4])^(alpha-1)-(1/(1+psi1))*x[4]^(1+psi1)-Cv/x[6]
    # h2:
    residuals[2]=alpha*1*ztilda(z,zlow,1-x[1])*(z*x[3]*x[4])^(alpha-1)-x[4]^psi1
    # q2:
    residuals[3]=x[6]-(mu*((1-x[3])./x[3]).^gamma).^(1/(1-gamma))
    
    # First period FOC:
    
    # n1:
    residuals[4]=alpha*epsilon*z*x[2]*(z*x[1]*x[2])^(alpha-1)-(1/(1+psi1))*x[2]^(1+psi1)-Cv/x[5]
    
    # h1:
    residuals[5]=alpha*epsilon*z*(z*x[1]*x[2])^(alpha-1)-x[2]^psi1
    
    # q1:
    residuals[6]=x[5]-(mu*((1-x[1])./x[1]).^gamma).^(1/(1-gamma))
    
end


Out[135]:
twoperiodeqflexible! (generic function with 1 method)

In [136]:
h1flexible=Array{Float64,1}(0)
h2flexible=Array{Float64,1}(0)
n1flexible=Array{Float64,1}(0)
n2flexible=Array{Float64,1}(0)

q1flexible=Array{Float64,1}(0)
q2flexible=Array{Float64,1}(0)

    # x[1]=n1; x[2]=h1; x[3]=n2z; x[4]=h2z; x[5]=n2zlow; x[6]=h2zlow; x[7]=q1; x[8]=q2z; x[9]=q2zlow

for eps=.5:.01:1
    epsilon=eps;
    r=nlsolve(twoperiodeqflexible!, [ .1; 1; .01; 1; .05; .05;])
    push!(n1flexible,r.zero[1])
    push!(h1flexible,r.zero[2])
    push!(n2flexible,r.zero[3])
    push!(h2flexible,r.zero[4])
    push!(q1flexible,r.zero[5])
    push!(q2flexible,r.zero[6])
end

In [137]:
hcat(n1flexible,n2flexible)


Out[137]:
51×2 Array{Float64,2}:
 0.451728  0.472196
 0.452811  0.472237
 0.453873  0.472277
 0.454916  0.472316
 0.455939  0.472355
 0.456944  0.472393
 0.457932  0.472431
 0.458903  0.472467
 0.459857  0.472503
 0.460795  0.472539
 0.461717  0.472574
 0.462625  0.472608
 0.463519  0.472641
 ⋮                 
 0.483453  0.473387
 0.484071  0.47341 
 0.484682  0.473433
 0.485287  0.473456
 0.485885  0.473478
 0.486477  0.4735  
 0.487063  0.473522
 0.487642  0.473543
 0.488216  0.473564
 0.488784  0.473585
 0.489346  0.473606
 0.489903  0.473627

In [138]:
hcat(h1flexible,h2flexible)


Out[138]:
51×2 Array{Float64,2}:
 0.519474  0.65524 
 0.525914  0.655545
 0.532308  0.655843
 0.538656  0.656137
 0.54496   0.656424
 0.55122   0.656707
 0.557439  0.656985
 0.563616  0.657257
 0.569753  0.657525
 0.575851  0.657789
 0.581911  0.658048
 0.587933  0.658303
 0.593919  0.658553
 ⋮                 
 0.744122  0.664133
 0.749332  0.664306
 0.75452   0.664476
 0.759687  0.664645
 0.764834  0.664812
 0.76996   0.664977
 0.775066  0.665141
 0.780152  0.665302
 0.785219  0.665462
 0.790266  0.665621
 0.795294  0.665778
 0.800304  0.665933

In [139]:
function twoperiodeqfixed!(x,residuals)
    # x[1]=n1; x[2]=n2; x[3]=q1; x[4]=q2; 
    
    # Last period FOC:
    
    # n2:
    residuals[1]=alpha*1*ztilda(z,zlow,1-x[1])*hbar2*(z*x[2]*hbar2)^(alpha-1)-(1/(1+psi1))*hbar2^(1+psi1)-Cv/x[4]

    # q2:
    residuals[2]=x[4]-(mu*((1-x[2])./x[2]).^gamma).^(1/(1-gamma))
    
    # First period FOC:
    
    # n1:
    residuals[3]=alpha*epsilon*z*hbar*(z*x[1]*hbar)^(alpha-1)-(1/(1+psi1))*hbar^(1+psi1)-Cv/x[3]
    
    # q1:
    residuals[4]=x[3]-(mu*((1-x[1])./x[1]).^gamma).^(1/(1-gamma))
    
end


Out[139]:
twoperiodeqfixed! (generic function with 1 method)

In [ ]:


In [140]:
n1fixed=Array{Float64,1}(0)
n2fixed=Array{Float64,1}(0)

q1fixed=Array{Float64,1}(0)
q2fixed=Array{Float64,1}(0)

    # x[1]=n1; x[2]=n2; x[3]=q1; x[4]=q2;
hbar=h1flexible[end]
hbar2=h2flexible[end]

for eps=.5:.01:1
    epsilon=eps;
    r=nlsolve(twoperiodeqfixed!, [ .1; .1; .05; .05;])
    push!(n1fixed,r.zero[1])
    push!(n2fixed,r.zero[2])
    push!(q1fixed,r.zero[3])
    push!(q2fixed,r.zero[4])
end

In [141]:
hcat(n1fixed,n2fixed)


Out[141]:
51×2 Array{Float64,2}:
 0.396415  0.468291
 0.40153   0.4686  
 0.406236  0.468883
 0.410583  0.469142
 0.414615  0.469381
 0.418368  0.469602
 0.421874  0.469808
 0.42516   0.470001
 0.42825   0.47018 
 0.431163  0.470349
 0.433917  0.470509
 0.436528  0.470659
 0.439007  0.470801
 ⋮                 
 0.480024  0.473092
 0.481018  0.473146
 0.48199   0.473199
 0.482942  0.473251
 0.483874  0.473301
 0.484787  0.473351
 0.485681  0.473399
 0.486558  0.473446
 0.487418  0.473493
 0.488262  0.473538
 0.48909   0.473583
 0.489903  0.473627

In [142]:
hcat(n1flexible,n2flexible)


Out[142]:
51×2 Array{Float64,2}:
 0.451728  0.472196
 0.452811  0.472237
 0.453873  0.472277
 0.454916  0.472316
 0.455939  0.472355
 0.456944  0.472393
 0.457932  0.472431
 0.458903  0.472467
 0.459857  0.472503
 0.460795  0.472539
 0.461717  0.472574
 0.462625  0.472608
 0.463519  0.472641
 ⋮                 
 0.483453  0.473387
 0.484071  0.47341 
 0.484682  0.473433
 0.485287  0.473456
 0.485885  0.473478
 0.486477  0.4735  
 0.487063  0.473522
 0.487642  0.473543
 0.488216  0.473564
 0.488784  0.473585
 0.489346  0.473606
 0.489903  0.473627

In [143]:
fig, ax = subplots()
ax[:plot](.5:.01:1,n1flexible./n1flexible[end], linewidth=2, color="blue", alpha=0.9, label="Period 1, flexible hours")
ax[:plot](.5:.01:1,n1fixed./n1fixed[end], linewidth=2, color="red", alpha=0.9, label="Period 1, fixed hours")

ax[:legend](loc="lower right")
ax[:set_title]("Employment in period 1")


Out[143]:
PyObject <matplotlib.text.Text object at 0x15bb56a10>

In [144]:
fig, ax = subplots()
ax[:plot](.5:.01:1,n2flexible./n2flexible[end], linewidth=2, color="blue", alpha=0.9, label="Period 2, flexible hours")
ax[:plot](.5:.01:1,n2fixed./n2fixed[end], linewidth=2, color="red", alpha=0.9, label="Period 2, fixed hours")

ax[:legend](loc="lower right")
ax[:set_title]("Employment in period 2")


Out[144]:
PyObject <matplotlib.text.Text object at 0x15bc13150>

In [150]:
output2flexible=(ztilda(z,zlow,1-n1flexible).*n2flexible.*h2flexible).^alpha;
output1flexible=collect(.5:.01:1).*(z*n1flexible.*h1flexible).^alpha;
output2fixed=(ztilda(z,zlow,1-n1fixed).*n2fixed.*hbar2).^alpha;
output1fixed=collect(.5:.01:1).*(z*n1fixed.*hbar).^alpha;

In [151]:
fig, ax = subplots()
ax[:plot](.5:.01:1,output1flexible./output1flexible[end], linewidth=2, color="blue", alpha=0.9, label="Period 1, flexible hours")
ax[:plot](.5:.01:1,output1fixed./output1fixed[end], linewidth=2, color="red", alpha=0.9, label="Period 1, fixed hours")

ax[:legend](loc="lower right")
ax[:set_title]("Output in period 1")


Out[151]:
PyObject <matplotlib.text.Text object at 0x15bddd690>

In [152]:
fig, ax = subplots()
ax[:plot](.5:.01:1,output2flexible./output2flexible[end], linewidth=2, color="blue", alpha=0.9, label="Period 2, flexible hours")
ax[:plot](.5:.01:1,output2fixed./output2fixed[end], linewidth=2, color="red", alpha=0.9, label="Period 2, fixed hours")

ax[:legend](loc="lower right")
ax[:set_title]("Output in period 2")


Out[152]:
PyObject <matplotlib.text.Text object at 0x15c167390>

In [156]:
fig, ax = subplots()
ax[:plot](.5:.01:1,h2flexible./h2flexible[end], linewidth=2, color="blue", alpha=0.9, label="Period 2, flexible hours")
ax[:plot](.5:.01:1,ones(length(.5:.01:1)), linewidth=2, color="red", alpha=0.9, label="Period 2, fixed hours")

ax[:legend](loc="lower right")
ax[:set_title]("Hours in period 2")


Out[156]:
PyObject <matplotlib.text.Text object at 0x15c43bd10>

In [155]:
fig, ax = subplots()
ax[:plot](.5:.01:1,h1flexible./h1flexible[end], linewidth=2, color="blue", alpha=0.9, label="Period 1, flexible hours")
ax[:plot](.5:.01:1,ones(length(.5:.01:1)), linewidth=2, color="red", alpha=0.9, label="Period 1, fixed hours")

ax[:legend](loc="lower right")
ax[:set_title]("Hours in period 1")


Out[155]:
PyObject <matplotlib.text.Text object at 0x156741890>

In [ ]: