chapter6.jl

# Chapter 6

#import Pkg

#Pkg.add("Plots")

#Pkg.add("Optim")

#Pkg.add("StatsKit")

#Pkg.add("GLM")

#Pkg.add("DataFrames")

#Pkg.add("StatsPlots")

#Pkg.add("StatFiles")

#Pkg.add("NLSolversBase")

#Pkg.add("Distributions")

using Random

using Statistics

using Plots

using Optim

using NLSolversBase

using StatsKit

using CSV

using GLM

using DataFrames

using LinearAlgebra

using StatsPlots

using StatFiles

using Distributions

# 6.1

# 6.2

# 6.2.1

# 6.2.2

rng = MersenneTwister(123457689)

N = 500

a = fill(2, N)

b = -3

x = rand!(rng, zeros(N))

υ = randn!(rng, zeros(N))

y_star = a + b.*x + υ

y = Vector{Union{Missing, Float64}}(missing,N)

for i = 1:N

if y_star[i] > 0

y[i] = y_star[i]

else

y[i] = 0

end

end

data1 = DataFrame(y = y, x = x)

lm1 = lm(@formula(y ~ x), data1)

coef(lm1)[2]

scatter(x, y)

scatter!(x, a + b.*x)

scatter!(x, coef(lm1)[1] .+ coef(lm1)[2].*x)

coef(lm(@formula(y ~ x), data1[data1.x .< 0.6,1:end]))[2]

length(data1.y[data1.x .< 0.6])

data2 = DataFrame(y = y .> 0, x = x)

coef(glm(@formula(y ~ x), data2, Binomial(), ProbitLink()))

# 6.2.3

histogram(y)

# 6.2.4

# 6.2.5

logϕ = function(z)

-z.^2 - √π

end

logΦ = function(v)

logcdf(Normal(), v)

end

Tobit = function(β, σ, y, X)

is0 = y .== 0

z = (y - X*β)./σ

return is0.*logΦ.(z) + (1 .- is0).*(logϕ.(z) .- log(σ))

end

Tobit_wrap = function(par)

σ = exp(par[1])

β = par[2:end]

return -sum(Tobit(β, σ, y, X))

end

global y = data1.y

global X = hcat(ones(N), data1.x)

init = [0; coef(lm1)]

a1 = optimize(Tobit_wrap, init)

exp(a1.minimizer[1])

a1.minimizer[2:end]

# 6.3

# 6.3.1

path = "/Volumes/My WD/Book"

x = DataFrame(CSV.read(string(path,"/NLSY97_min.csv")))

wage = Vector{Union{Missing,Float64}}(missing, 8984)

for i = 1:8984

if x[i,9] > 0 && x[i,8] > 0

wage[i] = x[i,8]/x[i,9]

end

end

# topcode wage

for i = 1:8984

if !ismissing(wage[i]) && wage[i] > quantile(skipmissing(wage), 0.9)

wage[i] = missing

end

end

lwage = Vector{Union{Missing,Float64}}(missing, 8984)

ed = Vector{Union{Missing,Float64}}(missing, 8984)

exper = Vector{Union{Missing,Float64}}(missing, 8984)

exper2 = Vector{Union{Missing,Float64}}(missing, 8984)

for i = 1:8984

if !ismissing(wage[i])

if wage[i] > 7.25

lwage[i] = log(wage[i])

else

lwage[i] = 0

end

if x.CV_HGC_EVER_EDT[i] > 0 && x.CV_HGC_EVER_EDT[i] < 25

ed[i] = x.CV_HGC_EVER_EDT[i]

end

exper[i] = 2010 - x[i,4] - ed[i] - 6

exper2[i] = (exper[i]^2)/100

end

println(i)

end

female = x.KEY!SEX .== 2

black = x.KEY!RACE_ETHNICITY .== 1

x1 = DataFrame(lwage=lwage, wage=wage, ed=ed, exper=exper, exper2=exper2,female=female,black=black, lwage_bin=lwage_bin)

x1 = dropmissing(x1)

x1.lwage_bin = x1.lwage .> 0

histogram(x1.wage[x1.wage .< 30])

vline!(fill(7.25,3972))

xlabel!("wage rate")

# 6.3.2

lm1 = lm(@formula(lwage ~ ed + exper + exper2 + female + black), x1)

glm1 = glm(@formula(lwage_bin ~ ed + exper + exper2 + female + black), x1, Binomial(), ProbitLink())

global y = x1.lwage

global X = hcat(ones(length(y)), x1.ed, x1.exper, x1.exper2, x1.female, x1.black)

init = [0; coef(glm1)]

a = optimize(Tobit_wrap, init)

exp(a.minimizer[1])

a.minimizer[2:end]

# 6.4

# 6.4.1

# 6.4.2

Random.seed!(123456789)

rng = MersenneTwister(123456789)

N = 100

a = fill(6, N)

b = -3

c = fill(4, N)

d = -5

x = rand!(rng, zeros(N))

z = rand!(rng, zeros(N))

μ = [0; 0]

ρ = 0.5

Σ = hcat([1; ρ], [ρ; 1])

Υ = rand(MvNormal(μ, Σ), N)'

y = Vector{Union{Missing, Float64}}(missing, N)

for i = 1:N

if c[i] + d*z[i] + Υ[i,1] > 0

y[i] = a[i] + b*x[i] + Υ[i,2]

else

y[i] = 0

end

end

# OLS Model

data1 = DataFrame(y = y, x = x)

lm1 = lm(@formula(y ~ x), data1)

# OLS Model with subset of data1

lm2 = lm(@formula(y ~ x), data1[z .< 0.6,1:end])

# Probit

data1.y1 = y .> 0

data1.z = z

glm1 = glm(@formula(y1 ~ z), data1, Binomial(), ProbitLink())

# 6.4.3

# 6.4.4

# 6.4.5

Heckman = function(β, γ, ρ, y, X, Z)

is0 = y .== 0

Zγ_adj = (Z*γ + ρ.*(y - X*β))./(√(1 - ρ^2))

return is0.*logΦ.(-Z*γ) + (1 .- is0).*(logΦ.(Zγ_adj) + logϕ.(y - X*β))

end

Heckman_wrap = function(par)

ρ = exp(par[1])/(1 + exp(par[2]))

J = size(X)[2]

β = par[2:(J+1)]

γ = par[(J+2):end]

return -sum(Heckman(β, γ, ρ, y, X, Z))

end

global y = data1.y

global X = hcat(ones(length(y)), data1.x)

global Z = hcat(ones(length(y)), data1.z)

init = [0; coef(lm1); coef(glm1)]

a1 = optimize(Heckman_wrap, init)

exp(a1.minimizer[1])/(1 + exp(a1.minimizer[1]))

a1.minimizer[2:3]

a1.minimizer[4:5]

# 6.5

# 6.5.1

path = "/Volumes/My WD/Book"

x = DataFrame(CSV.read(string(path,"/NLSY97_gender_book.csv")))

wage = Vector{Union{Missing,Float64}}(missing, 8984)

fulltime = Vector{Union{Missing,Float64}}(missing, 8984)

for i = 1:8984

if x[i,15] == "NA" || x[i,14] == "NA"

wage[i] = missing

elseif parse(Float64,x[i,15]) > 0

wage[i] = parse(Float64,x[i,14])/parse(Float64,x[i,15])

end

if x[i,15] == "NA"

fulltime[i] = missing

else

if parse(Float64,x[i,15]) > 1750

fulltime[i] = 1

else

fulltime[i] = 0

end

end

#println(i)

end

lwage = Vector{Union{Missing,Float64}}(missing, 8984)

lftwage = Vector{Union{Missing,Float64}}(missing, 8984)

age = Vector{Union{Missing,Float64}}(missing, 8984)

age2 = Vector{Union{Missing,Float64}}(missing, 8984)

college = Vector{Union{Missing,Float64}}(missing, 8984)

msa = Vector{Union{Missing,Float64}}(missing, 8984)

children = Vector{Union{Missing,Float64}}(missing, 8984)

for i = 1:8984

if !ismissing(wage[i])

if wage[i] > 1

lwage[i] = log(wage[i])

else

lwage[i] = 0

end

if lwage[i] > 0 && fulltime[i] == 1

lftwage[i] = lwage[i]

else

lftwage[i] = 0

end

if x.CV_HIGHEST_DEGREE_0708_2007[i] != "NA"

if parse(Int64, x.CV_HIGHEST_DEGREE_0708_2007[i]) > 2

college[i] = 1

else

college[i] = 0

end

end

msa[i] = x.CV_MSA_2007[i] ∈ ["2", "3", "4"]

if x.CV_BIO_CHILD_HH_2007[i] != "NA"

if parse(Int64,x.CV_BIO_CHILD_HH_2007[i]) > 0

children[i] = 1

else

children[i] = 0

end

end

end

println(i)

end

female = x.KEY_SEX_1997 .== 2

black = x.KEY_RACE_ETHNICITY_1997 .== 1

south = x.CV_CENSUS_REGION_2007 .== "3"

urban = x[1:end,13] .== "1"

married = x.CV_MARSTAT_COLLAPSED_2007 .== "2"

age = 2007 .- x.KEY_BDATE_Y_1997

age2 = age.^2

x1 = DataFrame(lftwage=lftwage,female=female, black=black, age=age, age2=age2, msa=msa, urban=urban, south=south, college=college, married=married, children=children)

x1 = dropmissing(x1)

x1.lftwage_bin = x1.lftwage .> 0

x1_f = x1[x1.female .== 1,1:end]

x1_m = x1[x1.female .== 0,1:end]

# Note that these samples are smaller than for the R code.

# not clear why there is a difference.

lm1 = lm(@formula(lftwage ~ age + age2 + black + college + south + msa), x1_m)

lm2 = lm(@formula(lftwage ~ age + age2 + black + college + south + msa), x1_f)

lm3 = lm(@formula(lftwage ~ female + age + age2 + black + college + south + msa), x1)

plot(kde(x1_m.lftwage[x1_m.lftwage_bin]))

plot!(kde(x1_f.lftwage[x1_f.lftwage_bin]))

# 6.5.2

glm1 = glm(@formula(lftwage_bin ~ college), x1_f, Binomial(), ProbitLink())

glm2 = glm(@formula(lftwage_bin ~ college + south), x1_f, Binomial(), ProbitLink())

glm3 = glm(@formula(lftwage_bin ~ college + south + married + children), x1_f, Binomial(), ProbitLink())

# 6.5.3

x1.lftwage_norm = (x1.lftwage .- mean(x1.lftwage))./std(x1.lftwage)

x1.lftwage_norm3 = x1.lftwage_norm

x1.lftwage_norm3[x1.lftwage .== 0] = zeros(sum(x1.lftwage .== 0))

global y = x1.lftwage_norm3[x1.female]

global X = hcat(ones(length(y)), x1.age[x1.female], x1.age2[x1.female], x1.black[x1.female], x1.college[x1.female], x1.south[x1.female], x1.msa[x1.female])

global Z = hcat(ones(length(y)), x1.college[x1.female], x1.south[x1.female], x1.married[x1.female], x1.children[x1.female])

init = [0; coef(lm2); coef(glm3)]

a = optimize(Heckman_wrap, init, show_trace=true, iterations=10000)

β̂ = a.minimizer[2:8]

yₐ = X*β̂ + randn!(rng, zeros(length(y)))

global y = x1.lftwage_norm3[x1.female .== 0]

global X = hcat(ones(length(y)), x1.age[x1.female .== 0], x1.age2[x1.female .== 0], x1.black[x1.female .== 0], x1.college[x1.female .== 0], x1.south[x1.female .== 0], x1.msa[x1.female .== 0])

y_am = X*β̂ + randn!(rng, zeros(length(y)))

plot(kde(yₐ))

plot!(kde(y_am))

plot!(kde(x1.lftwage_norm[(x1.female).*(x1.lftwage .> 0)]))

plot!(kde(x1.lftwage_norm[(x1.female .== 0).*(x1.lftwage .> 0)]))

# 6.6.1

path = "/Volumes/My WD/Book"

x = CSV.read(string(path,"/nls.csv"))

lwage76_temp = Vector{Union{Missing, Float64}}(missing,length(x.lwage76))

wage76_temp = Vector{Union{Missing, Float64}}(missing,length(x.wage76))

for i in 1:length(lwage76_temp)

if length(x.lwage76[i]) > 1

lwage76_temp[i] = parse(Float64, x.lwage76[i])

wage76_temp[i] = parse(Float64, x.wage76[i])

end

end

x.lwage76 = lwage76_temp

x.wage76 = wage76_temp

x1 = DataFrame(x)

x1 = dropmissing(x1)

x1.lwage76_norm = (x1.lwage76 .- mean(x1.lwage76))./std(x1.lwage76)

x1.exp1 = x1.age76 - x1.ed76 - fill(6,3010)

x1.exp2 = (x1.exp1.^2)./100

x1.college = x1.ed76 .> 12

# 6.6.2

lm_c = lm(@formula(lwage76_norm ~ exp1 + exp2 + black + reg76r + smsa76r), x1[x1.college,1:end])

lm_nc = lm(@formula(lwage76_norm ~ exp1 + exp2 + black + reg76r + smsa76r), x1[x1.college .== 0,1:end])

# 6.6.3

# Probit Estimate

glm1 = glm(@formula(college ~ nearc4), x1, Binomial(), ProbitLink())

glm2 = glm(@formula(college ~ nearc4 + momdad14), x1, Binomial(), ProbitLink())

glm3 = glm(@formula(college ~ nearc4 + momdad14 + black + smsa66r), x1, Binomial(), ProbitLink())

glm4 = glm(@formula(college ~ nearc4 + momdad14 + black + smsa66r + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + reg669), x1, Binomial(), ProbitLink())

# 6.6.4

# College

global y = x1.lwage76_norm[x1.college]

global X = hcat(ones(length(y)), x1.exp1[x1.college], x1.exp2[x1.college], x1.black[x1.college], x1.reg76r[x1.college], x1.smsa76r[x1.college])

global Z = hcat(ones(length(y)), x1.nearc4[x1.college], x1.momdad14[x1.college], x1.black[x1.college], x1.smsa66r[x1.college], x1.reg662[x1.college], x1.reg663[x1.college], x1.reg664[x1.college], x1.reg665[x1.college], x1.reg666[x1.college], x1.reg667[x1.college], x1.reg668[x1.college], x1.reg669[x1.college])

init = [0; coef(lm_c); coef(glm4)]

a_c = optimize(Heckman_wrap, init, show_trace=true, iterations=10000)

# No college

global y = -x1.lwage76_norm[x1.college .== 0]

global X = hcat(ones(length(y)), x1.exp1[x1.college .== 0], x1.exp2[x1.college .== 0], x1.black[x1.college .== 0], x1.reg76r[x1.college .== 0], x1.smsa76r[x1.college .== 0])

global Z = hcat(ones(length(y)), x1.nearc4[x1.college .== 0], x1.momdad14[x1.college .== 0], x1.black[x1.college .== 0], x1.smsa66r[x1.college .== 0], x1.reg662[x1.college .== 0], x1.reg663[x1.college .== 0], x1.reg664[x1.college .== 0], x1.reg665[x1.college .== 0], x1.reg666[x1.college .== 0], x1.reg667[x1.college .== 0], x1.reg668[x1.college .== 0], x1.reg669[x1.college .== 0])

init = [0; coef(lm_nc); -coef(glm4)]

a_nc = optimize(Heckman_wrap, init, show_trace=true, iterations=10000)

# Predicted wages

μ = [0; 0]

ρ̂ = exp(a_c.minimizer[1])/(1 + exp(a_c.minimizer[1]))

Σ̂ = [[1; ρ̂], [ρ̂; 1]]

Υ = rand(MvNormal(μ, Σ̂), size(x1)[1])'

β̂_c = a_c.minimizer[2:(length(coef(lm_c))+1)]

β̂_nc = a_nc.minimizer[2:(length(coef(lm_nc))+1)]

global X = hcat(ones(size(x1)[1]), x1.exp1, x1.exp2, x1.black, x1.reg76r, x1.smsa76r)

x1.coll_y = X*β̂_c + Υ[1:end,1]

x1.nocoll_y = -X*β̂_nc + Υ[1:end,2]

scatter(x1.nocoll_y, x1.coll_y)

scatter!(x1.nocoll_y, x1.nocoll_y)

xlabel!("wage (no-college)")

ylabel!("wage (college)")

# 6.6.5

mean(x1.coll_y - x1.nocoll_y)/4