chapter10.jl

# Chapter 10

#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

# 10.1

# 10.2

# 10.2.1

# 10.2.2

rng = MersenneTwister(123456789)

N = 1000

T = 2

a = fill(-2, N)

b = 3

d = -4

υ = [5; 6]

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

x₂ = x₁ + rand!(rng, zeros(N)) # change in X over time

x = hcat(x₁, x₂)'

e = randn!(rng, zeros(T, N))

y = hcat(a, a)' + b.*x .+ d.*υ + e

# 10.2.3

diffy = y[2,1:end] - y[1,1:end]

diffx = x[2,1:end] - x[1,1:end]

data1 = DataFrame(diffy=diffy, diffx=diffx)

lm1 = lm(@formula(diffy ~ diffx), data1)

# 10.3

# 10.3.1

# 10.3.2

did = function(y, treat)

y₁ = y[1,1:end]

y₂ = y[2,1:end]

did = Array{Union{Missing, String, Float64}}(missing, 4, 4)

did[1,2:4] = ["Period 2" "Period 1" "Diff"]

did[2:4,1] = ["Treated" "Not Treated" "Diff"]

did[2,2] = mean(y₂[treat .== 1])

did[3,2] = mean(y₁[treat .== 1])

did[2,3] = mean(y₂[treat .== 0])

did[3,3] = mean(y₁[treat .== 0])

did[4,2] = did[2,2] - did[3,2]

did[4,3] = did[2,3] - did[3,3]

did[2,4] = did[2,2] - did[2,3]

did[3,4] = did[3,2] - did[3,3]

did[4,4] = did[4,2] - did[4,3]

return did

end

rng = MersenneTwister(123456789)

treat = rand!(rng, zeros(N)) .< 0.5

x₂ = x₁ + treat

# this time the change is 1 for the treated group

x = hcat(x₁, x₂)'

y = hcat(a, a)' + b.*x .+ d.*υ + e

did1 = did(y, treat)

# 10.4

# 10.4.1

path = "/Volumes/My WD/Book"

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

EMPFT = Vector{Union{Missing, Float64}}(missing, 410)

EMPPT = Vector{Union{Missing, Float64}}(missing, 410)

EMPFT2 = Vector{Union{Missing, Float64}}(missing, 410)

EMPPT2 = Vector{Union{Missing, Float64}}(missing, 410)

for i = 1:410

if !(x.EMPFT[i] == "." || x.EMPPT[i] == "." || x.EMPFT2[i] == "." || x.EMPPT2[i] == ".")

EMPFT[i] = parse(Float64, x.EMPFT[i])

EMPPT[i] = parse(Float64, x.EMPPT[i])

EMPFT2[i] = parse(Float64, x.EMPFT2[i])

EMPPT2[i] = parse(Float64, x.EMPPT2[i])

end

x1 = DataFrame(SHEET = x.SHEET, FT = EMPFT, PT = EMPPT, STATE = x.STATE, TIME = ones(410))

x1.SHEET[x1.SHEET .== 407] = [4071; 4072]

# there is an issue with the labeling in the data

# two firms with the same label

x2 = DataFrame(SHEET = x.SHEET, FT = EMPFT2, PT = EMPPT2, STATE = x.STATE, TIME = ones(410) .+ 1)

x2.SHEET[x2.SHEET .== 407] = [4071; 4072]

x3 = vcat(x1, x2)

x4 = dropmissing!(x3)

wage = []

for i = 1:410

if !(x.WAGE_ST[i] == ".")

wage = [wage; parse(Float64, x.WAGE_ST[i])]

end

histogram(wage, legend = :none)

vline!(fill(4.25, 410))

vline!(fill(5.05, 410))

xlabel!("Wage")

ylabel!("Frequency")

# 10.4.2

y = hcat(x4.FT[x4.TIME .== 1] + x4.PT[x4.TIME .== 1], x4.FT[x4.TIME .== 2] + x4.PT[x4.TIME .== 2])

treat = x4.STATE[x4.TIME .== 1] .== 1

did2 = did(y', treat) # note that y is y'

# 10.5

# 10.5.1

# 10.5.2

# 10.5.3

# 10.5.4

# 10.5.5

fe = function(Y; X = false, cross = true)

T = size(Y)[1]

N = size(Y)[2]

y = Y[1:end,1]

t = 1:T

for i = 2:N

y = vcat(y, Y[1:end,i])

t = vcat(t, 1:T)

end

# set up for different cases

if cross

c = fill(1, T)

for i = 2:N

c = vcat(c, fill(i, T))

end

if length(X) > 1

treat = X[1:end,1]

for i = 2:N

treat = vcat(treat, X[1:end,i])

end

# estimator

if cross && length(X) > 1

data1 = DataFrame(y = y, treat = treat, t = string.(t), c = string.(c))

lm1 = lm(@formula(y ~ treat + t + c), data1)

elseif length(X) > 1

data1 = DataFrame(y = y, treat = treat, t = string.(t))

lm1 = lm(@formula(y ~ treat + t), data1)

else

data1 = DataFrame(y = y, t = string.(t))

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

end

return lm1

end

rng = MersenneTwister(123456789)

N = 100

T = 10

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

γ = rand!(rng, zeros(T))

β = 3

ϵ = randn!(rng, zeros(T, N))

treat = rand!(rng, zeros(N)) .< 0.5

y = fill(α[1], T) + γ + ϵ[1:end,1]

for i = 2:N

y = hcat(y, fill(α[i], T) + γ + ϵ[1:end,i])

end

y[1,1:end] = y[1,1:end] + β.*treat

treat1 = zeros(T, N)

treat1[1,1:end] = treat

# standard estimator

lm1 = fe(y, X=treat1)

coef(lm1)[2]

# No individual fixed effect estimator

lm2 = fe(y, X=treat1, cross=false)

coef(lm2)[2]

# Adjusted estimator

y₀ = y[2:T,1:end]

α̂ = fill(mean(y[1:end,1]), T)

for i = 2:N

α̂ = hcat(α̂, fill(mean(y[1:end,i]), T))

end

y₂ = y - α̂

lm3 = fe(y₂, X=treat1, cross=false)

coef(lm3)[2]

# Two step estimator

lm4 = fe(y₀, cross=false) # adjust for time effects

α̂ = fill(mean(residuals(lm4)[1:(T-1)]), T)

for i = 2:N

α̂ᵢ = fill(mean(residuals(lm4)[(1 + (i-1)*(T-1)):(i*(T-1))]), T)

α̂ = hcat(α̂, α̂ᵢ)

end

y₃ = y - α̂

lm5 = fe(y₃, X=treat1, cross=false)

coef(lm5)[2]

# 10.6

# 10.6.1

path = "ENTER PATH HERE"

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

x_names = CSV.read(string(path,"/NLSY97Panel_names.csv"), header=false)

year = 1997:2014

year1 = [string.(97:99); "00"; "01"; "02"; "03"; "04"; "05"; "06"; "07"; "08"; "09"; string.(10:14)]

W = Array{Union{Missing, Float64}}(missing, 18, 8984)

Y = Array{Union{Missing, Float64}}(missing, 18, 8984)

for i = 1:18

hrs_name = string("CVC_HOURS_WK_YR_ALL_",year1[i],"_XRND")

inc_name = string("YINC_1700_",year[i])

for j = 1:148

if x_names[1][j] == hrs_name

for k = 1:8984

if x[k,j] >= 0

Y[i,k] = x[k,j]

end

if x_names[1][j] == inc_name

for k = 1:8984

if x[k,j] >= 0

W[i,k] = x[k,j]

end

rate_07 = W[11,1:end]./Y[11,1:end]

treat = Vector{Union{Missing, Int64}}(missing, 8984)

for i = 1:8984

if !(ismissing(rate_07[i]))

if rate_07[i] < 7.26

treat[i] = 1

else

treat[i] = 0

end

if !ismissing(W[11,i])

if W[11,i] == 0

treat[i] = 1

end

if !ismissing(Y[11,i])

if Y[11,i] == 0

treat[i] = 1

end

# 10.6.2

N = 600 # reduce the size for computational reasons

T = 18

treat1 = Array{Union{Missing, Int64}}(missing, T, N)

treat1[1:10,1:end] .= 0

for i = 11:T

treat1[i,1:end] = treat[1:N]

end

lm1 = fe(Y[1:end,1:N], X=treat1)

coef(lm1)[2]

N = 8984 # full data set

y₀ = Y[1:10,1:N]

α̂ = Array{Union{Missing, Float64}}(missing, T, N)

for i = 1:N

α̂[1:end,i] = fill(mean(skipmissing(y₀[1:end,i])), T)

end

y₁ = Y - α̂

treat1 = Array{Union{Missing, Int64}}(missing, T, N)

treat1[1:10,1:end] .= 0

for i = 11:T

treat1[i,1:end] = treat[1:N]

end

lm2 = fe(y₁, X=treat1, cross=false)

coef(lm2)[2]

lm3 = fe(y₀, cross=false)

y₀_res = y₀[1,1:end] .- coef(lm3)[1]

y₀_res = hcat(y₀_res, y₀[2,1:end] .- coef(lm3)[3])

y₀_res = hcat(y₀_res, y₀[3,1:end] .- coef(lm3)[4])

y₀_res = hcat(y₀_res, y₀[4,1:end] .- coef(lm3)[5])

y₀_res = hcat(y₀_res, y₀[5,1:end] .- coef(lm3)[6])

y₀_res = hcat(y₀_res, y₀[6,1:end] .- coef(lm3)[7])

y₀_res = hcat(y₀_res, y₀[7,1:end] .- coef(lm3)[8])

y₀_res = hcat(y₀_res, y₀[8,1:end] .- coef(lm3)[9])

y₀_res = hcat(y₀_res, y₀[9,1:end] .- coef(lm3)[10])

y₀_res = hcat(y₀_res, y₀[10,1:end] .- coef(lm3)[2])

y₀_res = y₀_res'

α̂ = fill(mean(skipmissing(y₀_res[1:end,1])), T)

for i = 2:N

α̂ᵢ = fill(mean(skipmissing(y₀_res[1:end,i])), T)

α̂ = hcat(α̂, α̂ᵢ)

end

y₃ = Y - α̂

lm4 = fe(y₃, X=treat1, cross=false)

coef(lm4)[2]

W_nm = W[1:end, ismissing.(treat) .== 0]

Y_nm = Y[1:end, ismissing.(treat) .== 0]

treat_nm = treat[ismissing.(treat) .== 0]

wage_07_nm = W_nm[11,treat_nm .== 1]./Y_nm[11,treat_nm .== 1]

for i = 1:length(wage_07_nm)

if !(ismissing(wage_07_nm[i]))

if isnan(wage_07_nm[i]) || Y_nm[11,treat_nm .== 1][i] == 0

wage_07_nm[i] = missing

end

c = mean(skipmissing(wage_07_nm)) # 2007 wage rate

d = mean(skipmissing(W_nm[14,treat_nm .== 1])) # 2010 income

e = mean(skipmissing(Y_nm[14,treat_nm .== 1])) # 2010 hours

d - c*(e + 270) # actual less counterfactual

# Note we may not be accounting for missing the same way as the R code.

# Also I calculated the average of the wage rate for the treated

# not the average income divided by average hours for the treated.

# 10.7