In [325]:

    

using NLsolve
using PyPlot


    

In [326]:

    

# model parameters
Cv=0.2;
psi1=1.01;
alpha=0.5;
epsilon=1.0;
hbar=1.0;
gamma=.85;
mu=0.9;
z=1;
zlow=.5;


    

In [327]:

    

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


    

    
Out[327]:





q (generic function with 1 method)

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [328]:

    

function twoperiodeqflexiblehours!(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; x[9]=q2zlow
    
    # Last period FOC:
    
    # n2z:
    residuals[1]=alpha*1*z*x[4]*(z*x[3]*x[4])^(alpha-1)-(1/(1+psi1))*x[4]^(1+psi1)-Cv/x[8]
    # n2zlow:
    residuals[2]=alpha*1*zlow*x[6]*(zlow*x[5]*x[6])^(alpha-1)-(1/(1+psi1))*x[6]^(1+psi1)-Cv/x[9]
    # h2z:
    residuals[3]=alpha*1*z*(z*x[3]*x[4])^(alpha-1)-x[4]^psi1
    # h2zlow:
    residuals[4]=alpha*1*zlow*(zlow*x[5]*x[6])^(alpha-1)-x[6]^psi1
    # q2z:
    residuals[5]=x[8]-(mu*((x[1]-x[3])./x[3]).^gamma).^(1/(1-gamma))
    # q2zlow:
    residuals[6]=x[9]-(mu*((1-x[1]-x[5])./x[5]).^gamma).^(1/(1-gamma))
    
    # First period FOC:
    
    # n1:
    residuals[7]=alpha*epsilon*z*x[2]*(z*x[1]*x[2])^(alpha-1)-(1/(1+psi1))*x[2]^(1+psi1)-Cv/x[7]-(Cv/x[9])*(gamma/(1-gamma))*(x[5]/(1-x[1]-x[5]))+(Cv/x[8])*(gamma/(1-gamma))*(x[3]/(x[1]-x[3]))
    
    # h1:
    residuals[8]=alpha*epsilon*z*(z*x[1]*x[2])^(alpha-1)-x[2]^psi1
    
    # q1:
    residuals[9]=x[7]-(mu*((1-x[1])./x[1]).^gamma).^(1/(1-gamma))
    
end


    

    
Out[328]:





twoperiodeqflexiblehours! (generic function with 1 method)

In [329]:

    

h1flexible=Array{Float64,1}(0)
h2zflexible=Array{Float64,1}(0)
h2zlowflexible=Array{Float64,1}(0)

n1flexible=Array{Float64,1}(0)
n2zflexible=Array{Float64,1}(0)
n2zlowflexible=Array{Float64,1}(0)

q1flexible=Array{Float64,1}(0)
q2zflexible=Array{Float64,1}(0)
q2zlowflexible=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=.1:.1:1
    epsilon=eps;
    r=nlsolve(twoperiodeqflexiblehours!, [ .1; 1; .01; 1; .005; 1; .05; .05; .05])
    push!(n1flexible,r.zero[1])
    push!(h1flexible,r.zero[2])
    push!(n2zflexible,r.zero[3])
    push!(h2zflexible,r.zero[4])
    push!(n2zlowflexible,r.zero[5])
    push!(h2zlowflexible,r.zero[6])
    push!(q1flexible,r.zero[7])
    push!(q2zflexible,r.zero[8])
    push!(q2zlowflexible,r.zero[9])
end


    

In [330]:

    

hcat(h1flexible,h2zflexible,h2zlowflexible)


    

    
Out[330]:





10×3 Array{Float64,2}:
 0.16898   0.972129  0.819696
 0.267294  0.971701  0.820115
 0.349438  0.971202  0.820605
 0.422525  0.970655  0.821145
 0.489504  0.970073  0.821721
 0.551959  0.969467  0.822327
 0.610866  0.968842  0.822954
 0.666881  0.968203  0.823599
 0.72047   0.967554  0.824258
 0.771983  0.966899  0.824927

In [331]:

    

hcat(n1flexible,n2zflexible,n2zlowflexible)


    

    
Out[331]:





10×3 Array{Float64,2}:
 0.536882  0.272279  0.227865
 0.537639  0.272641  0.227514
 0.538522  0.273065  0.227104
 0.539492  0.27353   0.226654
 0.540525  0.274025  0.226174
 0.541607  0.274543  0.225671
 0.542724  0.275079  0.225152
 0.543869  0.275627  0.22462 
 0.545035  0.276185  0.224078
 0.546215  0.276751  0.223529

In [332]:

    

hcat(q1flexible,q2zflexible,q2zlowflexible)


    

    
Out[332]:





10×3 Array{Float64,2}:
 0.214397  0.421289  0.593557
 0.210733  0.421662  0.592948
 0.206532  0.422098  0.592236
 0.202013  0.422576  0.591454
 0.197306  0.423085  0.59062 
 0.192497  0.423617  0.589747
 0.187648  0.424167  0.588843
 0.182805  0.42473   0.587917
 0.178001  0.425302  0.586973
 0.173262  0.425881  0.586015

In [333]:

    

function twoperiodeqfixedhours!(x,residuals)
    # x[1]=n1; x[2]=n2z; x[3]=n2zlow; x[4]=q1; x[5]=q2z; x[6]=q2zlow
    
    # Last period FOC:
    
    # n2z:
    residuals[1]=alpha*1*z*hbarz*(z*x[2]*hbarz)^(alpha-1)-(1/(1+psi1))*hbarz^(1+psi1)-Cv/x[5]
    # n2zlow:
    residuals[2]=alpha*1*zlow*hbarzlow*(zlow*x[3]*hbarzlow)^(alpha-1)-(1/(1+psi1))*hbarzlow^(1+psi1)-Cv/x[6]
    # q2z:
    residuals[3]=x[5]-(mu*((x[1]-x[2])./x[2]).^gamma).^(1/(1-gamma))
    # q2zlow:
    residuals[4]=x[6]-(mu*((1-x[1]-x[3])./x[3]).^gamma).^(1/(1-gamma))
    
    # First period FOC:
    
    # n1:
    residuals[5]=alpha*epsilon*z*hbar*(z*x[1]*hbar)^(alpha-1)-(1/(1+psi1))*hbar^(1+psi1)-Cv/x[4]-(Cv/x[6])*(gamma/(1-gamma))*(x[3]/(1-x[1]-x[3]))+(Cv/x[5])*(gamma/(1-gamma))*(x[2]/(x[1]-x[2]))
    
    # q1:
    residuals[6]=x[4]-(mu*((1-x[1])./x[1]).^gamma).^(1/(1-gamma))
    
end


    

    
Out[333]:





twoperiodeqfixedhours! (generic function with 1 method)

In [334]:

    

hbar=h1flexible[end];
hbarz=h2zflexible[end];
hbarzlow=h2zlowflexible[end];
n1fixed=Array{Float64,1}(0)
n2zfixed=Array{Float64,1}(0)
n2zlowfixed=Array{Float64,1}(0)

q1fixed=Array{Float64,1}(0)
q2zfixed=Array{Float64,1}(0)
q2zlowfixed=Array{Float64,1}(0)

# x[1]=n1; x[2]=n2z; x[3]=n2zlow; x[4]=q1; x[5]=q2z; x[6]=q2zlow

for eps=.1:.1:1
    epsilon=eps;
    r=nlsolve(twoperiodeqfixedhours!, [ .1; .01; .1; .05; .05; .05])
    push!(n1fixed,r.zero[1])
    push!(n2zfixed,r.zero[2])
    push!(n2zlowfixed,r.zero[3])
    push!(q1fixed,r.zero[4])
    push!(q2zfixed,r.zero[5])
    push!(q2zlowfixed,r.zero[6])
end


    

In [335]:

    

hcat(n1fixed,n2zfixed,n2zlowfixed)


    

    
Out[335]:





10×3 Array{Float64,2}:
 0.529369  0.268881  0.231131
 0.531502  0.269879  0.230171
 0.533562  0.270842  0.229243
 0.535551  0.271772  0.228346
 0.537475  0.272671  0.227478
 0.539336  0.273541  0.226638
 0.541137  0.274382  0.225824
 0.542883  0.275196  0.225036
 0.544574  0.275986  0.224271
 0.546215  0.276751  0.223529

In [336]:

    

hcat(n1flexible,n2zflexible,n2zlowflexible)


    

    
Out[336]:





10×3 Array{Float64,2}:
 0.536882  0.272279  0.227865
 0.537639  0.272641  0.227514
 0.538522  0.273065  0.227104
 0.539492  0.27353   0.226654
 0.540525  0.274025  0.226174
 0.541607  0.274543  0.225671
 0.542724  0.275079  0.225152
 0.543869  0.275627  0.22462 
 0.545035  0.276185  0.224078
 0.546215  0.276751  0.223529

In [337]:

    

hcat(h1flexible,h2zflexible,h2zlowflexible)


    

    
Out[337]:





10×3 Array{Float64,2}:
 0.16898   0.972129  0.819696
 0.267294  0.971701  0.820115
 0.349438  0.971202  0.820605
 0.422525  0.970655  0.821145
 0.489504  0.970073  0.821721
 0.551959  0.969467  0.822327
 0.610866  0.968842  0.822954
 0.666881  0.968203  0.823599
 0.72047   0.967554  0.824258
 0.771983  0.966899  0.824927

In [338]:

    

output2flexible=(z*n2zflexible.*h2zflexible).^alpha+(zlow*n2zlowflexible.*h2zlowflexible).^alpha;
output2fixed=(z*n2zfixed.*hbarz).^alpha+(zlow*n2zlowfixed.*hbarzlow).^alpha;
output1flexible=collect(.1:.1:1).*(z*n1flexible.*h1flexible).^alpha;
output1fixed=collect(.1:.1:1).*(z*n1fixed.*hbar).^alpha;
hcat(output1flexible,output2flexible)


    

    
Out[338]:





10×2 Array{Float64,2}:
 0.0301201  0.820079
 0.0758176  0.82015 
 0.130139   0.820233
 0.190976   0.820323
 0.257191   0.820419
 0.328055   0.820518
 0.403052   0.82062 
 0.481794   0.820723
 0.563979   0.820827
 0.649361   0.820932

In [339]:

    

hcat(output1fixed,output2fixed)


    

    
Out[339]:





10×2 Array{Float64,2}:
 0.0639268  0.818644
 0.128111   0.818947
 0.192538   0.819236
 0.257196   0.819512
 0.322072   0.819775
 0.387155   0.820027
 0.452435   0.820268
 0.517901   0.820499
 0.583546   0.82072 
 0.649361   0.820932

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: