chapter11.jl

# Chapter 11

#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

# 11.1

# 11.2

# 11.2.1

# 11.2.2

# 11.2.3

rng = MersenneTwister(123456789)

N = 200 # individuals

T = 35 # time periods

K = 3 # factors

σ = 0.1 # noise

F = rand!(rng, zeros(T, K)).*3

Λ = exp.(randn!(rng, zeros(K, N)))

for i = 1:N

col_sum = sum(Λ[1:end,i])

Λ[1:end,i] = Λ[1:end,i]./col_sum

end

E = randn!(rng, zeros(T, N)).*σ

Ȳ = F*Λ

Ys = Ȳ + E

scatter(Ys[2,1:end], Ys[1,1:end])

scatter!(F[2,1:end], F[1,1:end])

xlabel!("Period 2")

ylabel!("Period 1")

# 11.2.4

# 11.2.5

# 11.2.6

tab_res = Array{Union{Missing, Float64}}(missing, 34, 2)

for i = 2:35

y = Ys[1:30,1]

X = hcat(ones(30), Ys[1:30,2:i])

ω̂ = inv(X'X)X'y

res = Ys[31:35,1] - hcat(ones(5),Ys[31:35,2:i])*ω̂

tab_res[i-1,1] = i

tab_res[i-1,2] = √(mean(res.^2))

end

plot(tab_res[1:end,1],tab_res[1:end,2], legend = :none)

ylabel!("Root Mean Squared Error")

xlabel!("Non-Treated Units Used")

y = Ys[1:30,1]

X = hcat(ones(30), Ys[1:30,2:7])

ω̂ = inv(X'X)X'y

ŷ_ols = hcat(ones(5), Ys[31:35,2:7])*ω̂

res_ols = √mean((Ys[31:35,1] - ŷ_ols).^2)

# Not sure R could actually solve this.

# 11.3

# 11.3.1

scatter(Ys[1,1:end], Ys[2,1:end])

scatter!(Ys[1,1:2],Ys[2,1:2])

# 11.3.2

sc = function(par; yt=false, ynt=false)

ω = exp.(par)./sum(exp.(par))

return mean((yt - ynt*ω).^2)

end

J = size(X[1:end,2:end])[2]

init = randn!(rng, zeros(J))

par = init

a = optimize(par -> sc(par, yt = Ys[1:30,1], ynt = X[1:30,2:end]), init, iterations = 1000000, show_trace = false)

ω̂ = exp.(a.minimizer)./sum(exp.(a.minimizer))

ŷ_sc = Ys[31:35,2:end]*ω̂

res_sc = √mean((Ys[31:35,1] - ŷ_sc).^2)

# 11.4

# 11.4.1

# 11.4.2

# 11.4.3

ols_reg = function(ω; yt = false, Ynt = false, θ = 0)

X = hcat(ones(length(yt)), Ynt)

return mean((yt - X*ω).^2) + θ*sum(abs.(ω))

end

y_r = Ys[1:25,1] # training data

y_t = Ys[26:30,1] # test data

Y_r = Ys[1:25,2:end]

Y_t = Ys[26:30,2:end]

lasso = function(y_r, y_t, Y_r, Y_t; K = 3)

θ_values = 0:(K-1)

J = size(Y_r)[2]

global ω_th = zeros(J+1)

global rmse_old = 10000 # some large number

for k = 1:K

init = zeros(J+1)

a_k = optimize(par -> ols_reg(par, yt = y_r, Ynt = Y_r, θ = θ_values[k]), init, iterations=1000000)

ω̂ = a_k.minimizer

X = hcat(ones(length(y_t)), Y_t)

global rmse_new = √mean((y_t - X*ω̂).^2)

if rmse_new < rmse_old

global ω_th = ω̂

global rmse_old = rmse_new

println(rmse_new)

end

return ω_th

end

ω̂ = lasso(y_r, y_t, Y_r, Y_t, K=5)

X = hcat(ones(5), Ys[31:35,2:end])

ŷ_lasso = X*ω̂

res_lasso = √mean((Ys[31:35,1] - ŷ_lasso).^2)

# 11.5

# 11.5.1

# 11.5.2

# 11.5.3

# 11.5.4

# 11.5.5

logϕ = function(z)

-0.5*z.^2 - √(2*π)

end

ϕ = function(z)

(1/√(2*π))*exp(-0.5*z^2)

end

cmf = function(par, Y, K; Reps=10)

ϵ = 1e-20

rng = MersenneTwister(123456789)

T = size(Y)[1]

N = size(Y)[2]

σ = exp(par[1])

F = zeros(T, K)

for k = 1:K

F[1:end,k] = par[(2 + (k-1)*T):(1 + k*T)]

end

F = exp.(F)

p = zeros(T, N)

for r = 1:Reps

L = rand!(rng, zeros(K,N))

for i = 1:N

colsum = sum(L[1:end,i])

L[1:end,i] = L[1:end,i]./colsum

end

p = p + ϕ.((Y - F*L)./σ)

end

p = p./Reps

return - (mean(log.(p .+ ϵ)) - log(σ))

end

init = [log(0.1); randn!(rng, zeros(15))]

a = optimize(par -> cmf(par, Ys[1:5,1:end], 3), init, show_trace=false, iterations=100000)

σ̂ = exp(a.minimizer[1])

T = 5

F̂ = zeros(T, K)

for k = 1:K

F̂[1:end,k] = a.minimizer[(2 + (k-1)*T):(1 + k*T)]

end

F̂ = exp.(F̂)

λ̂ = function(par, Yᵢ, F̂, σ̂)

K = size(F̂)[1]

L = exp.(par)./sum(exp.(par))

return - sum(logϕ.((Yᵢ - F̂*L)./σ̂))

end

rng = MersenneTwister(123456789)

Λ̂ = zeros(3, 200)

for i = 1:200

init = randn!(rng, zeros(K)).*20

aᵢ = optimize(par -> λ̂(par, Ys[1:5,i], F̂, σ̂), init)

Λ̂[1:end,i] = exp.(aᵢ.minimizer)./sum(exp.(aᵢ.minimizer))

#println(i)

end

scatter(Λ̂[1,1:end], Λ̂[2,1:end])

xlabel!("Factor 1")

ylabel!("Factor 2")

Λ̂ⱼ = Λ̂[1:end,2:200]

λ̂ᵢ = Λ̂[1:end,1]

Yⱼ = Ys[31:35,2:200]

ŷ_fm = Yⱼ*Λ̂ⱼ'*inv(Λ̂ⱼ*Λ̂ⱼ')*λ̂ᵢ

res_fm = √mean((Ys[31:35,1] - ŷ_fm).^2)

plot(Ys[31:35,1],Ys[31:35,1], label = :none)

scatter!(Ys[31:35,1], ŷ_ols, label = "OLS")

scatter!(Ys[31:35,1], ŷ_sc, label = "SC")

scatter!(Ys[31:35,1], ŷ_lasso, label = "LASSO")

scatter!(Ys[31:35,1], ŷ_fm, label = "FM", legend = :topleft)

# 11.6

# 11.6.1

path = "ENTER PATH"

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

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

Y2 = Y[1:end,1:N]

W2 = W[1:end,1:N]

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

lm1 = fe(Y2, X=treat1)

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

Y3[(ismissing.(Y2) .== 0).*(ismissing.(treat1) .== 0)] = residuals(lm1)

Y_t = Array{Union{Missing, Float64}}(missing, T, sum(skipmissing(treat)))

k = 1

for i = 1:N

if !ismissing(treat[i])

if treat[i] == 1

Y_t[1:end,k] = Y3[1:end,i]

k += 1

end

Y_t2 = Y_t[1:end,2]

for i = 2:N

if !ismissing(sum(Y_t[1:end,i]))

Y_t2 = hcat(Y_t2, Y_t[1:end,i])

end

Y_t2 = Y_t2[1:end,1:5]

Y_nt = Array{Union{Missing, Float64}}(missing, T, sum(skipmissing(1 .- treat)))

k = 1

for i = 1:N

if !ismissing(treat[i])

if treat[i] == 0

Y_nt[1:end,k] = Y3[1:end,i]

k += 1

end

Y_nt2 = Y_nt[1:end,2]

for i = 2:N

if !ismissing(sum(Y_nt[1:end,i]))

Y_nt2 = hcat(Y_nt2, Y_nt[1:end,i])

end

# 11.6.2

rng = MersenneTwister(123456789)

N_t = size(Y_t2)[2]

N_nt = size(Y_nt2)[2]

mat_treat = Array{Union{Missing, Float64}}(missing, 8, N_t)

for i = 1:N_t

init = randn!(rng, zeros(N_nt))

aᵢ = optimize(par -> sc(par, yt=Y_t2[1:10,i], ynt = Y_nt2[1:10,1:end]), init, show_trace=false, iterations=1000000)

ω̂ = exp.(aᵢ.minimizer)./sum(exp.(aᵢ.minimizer))

mat_treat[1:end,i] = Y_t2[11:18,i] - Y_nt2[11:18,1:end]*ω̂

println(i)

end

mat_treat2 = mat_treat .+ coef(lm1)[2]

# add back the fixed effects results

plot_res_sc = Array{Union{Missing, Float64}}(missing, 8, 3)

for i = 1:8

plot_res_sc[i,1:end] = quantile(mat_treat2[i,1:end], [0.2; 0.5; 0.8])

end

# 11.6.3

Y_tr2 = Y_t2[1:8,1:end]

Y_tt2 = Y_t2[9:10,1:end]

Y_ntr2 = Y_nt2[1:8,1:end]

Y_ntt2 = Y_nt2[9:10,1:end]

treat3 = Array{Union{Missing, Float64}}(missing, 8, 5)

for i = 1:5

ω̂ = lasso(Y_tr2[1:end,i], Y_tt2[1:end,i], Y_ntr2, Y_ntt2, K=10)

treat3[1:end,i] = Y_t2[11:18,i] - hcat(ones(8), Y_nt2[11:18,1:end])*ω̂

println(i)

end

treat3 = treat3 .+ coef(lm1)[2]

plot_res_ls = Array{Union{Missing, Float64}}(missing, 8, 3)

for i = 1:8

plot_res_ls[i,1:end] = quantile(treat3[i,1:end], [0.2; 0.5; 0.8])

end

plot(2009:2016, plot_res_sc[1:end,1], lw=3, color = :black, label = "SC")

plot!(2009:2016, plot_res_sc[1:end,2], lw=3, color = :black, label = :none)

plot!(2009:2016, plot_res_sc[1:end,3], lw=3, color = :black, label = :none)

plot!(2009:2016, plot_res_ls[1:end,1], lw=3, color = :black, style = :dash, label = "LASSO")

plot!(2009:2016, plot_res_ls[1:end,2], lw=3, color = :black, style = :dash, label = :none)

plot!(2009:2016, plot_res_ls[1:end,3], lw=3, color = :black, style = :dash, label = :none)

xlabel!("Year")

ylabel!("Change in Hours")

# 11.6.4

rng = MersenneTwister(123456789)

Y_t = Array{Union{Missing, Float64}}(missing, 18, sum(skipmissing(treat)))

W_t = Array{Union{Missing, Float64}}(missing, 18, sum(skipmissing(treat)))

Y_nt = Array{Union{Missing, Float64}}(missing, 18, sum(skipmissing(1 .- treat)))

k1 = 1

k2 = 1

for i = 1:8984

if !ismissing(treat[i])

if treat[i] == 1

Y_t[1:end,k1] = Y[1:end,i]

W_t[1:end,k1] = W[1:end,i]

k1 += 1

else

Y_nt[1:end,k2] = Y[1:end,i]

k2 += 1

end

Y_t2 = Y_t[1:end,2]

index_t_na = [2]

for i = 2:(size(Y_t)[2])

if !ismissing(sum(Y_t[1:end,i]))

Y_t2 = hcat(Y_t2, Y_t[1:end,i])

index_t_na = [index_t_na; i]

end

Y_t3 = zeros(3,1545)

for i = 1:1545

Y_t3[1,i] = mean(Y_t2[1:5,i])

Y_t3[2,i] = mean(Y_t2[6:10,i])

Y_t3[3,i] = mean(Y_t2[11:15,i])

end

Y_nt2 = Y_nt[1:end,2]

for i = 2:(size(Y_nt)[2])

if !ismissing(sum(Y_nt[1:end,i]))

Y_nt2 = hcat(Y_nt2, Y_nt[1:end,i])

end

Y_nt3 = zeros(3,2606)

for i = 1:2606

Y_nt3[1,i] = mean(Y_nt2[1:5,i])

Y_nt3[2,i] = mean(Y_nt2[6:10,i])

Y_nt3[3,i] = mean(Y_nt2[11:15,i])

end

Y_fm = hcat(Y_t3, Y_nt3)

N_t = 1545

N_nt = 2206

T = 2

K = 3

init = [log(400); randn!(rng, zeros(K*T)).*10]

a = optimize(par -> cmf(par, Y_fm[1:2,1:end], K, Reps=100), init, show_trace=true, iterations=100000)

σ̂ = exp(a.minimizer[1])

F̂ = zeros(T, K)

for k = 1:K

F̂[1:end,k] = a.minimizer[(2 + (k-1)*T):(1 + k*T)]

end

F̂ = exp.(F̂)

rng = MersenneTwister(123456789)

Λ̂ = zeros(3, 4151)

for i = 1:4151

init = randn!(rng, zeros(K)).*5

aᵢ = optimize(par -> λ̂(par, Y_fm[1:2,i], F̂, σ̂), init)

Λ̂[1:end,i] = exp.(aᵢ.minimizer)./sum(exp.(aᵢ.minimizer))

#println(i)

end

scatter(Λ̂[1,1:end], Λ̂[2,1:end])

xlabel!("Factor 1")

ylabel!("Factor 2")

L̂_t = Λ̂[1:end,1:N_t]

L̂_nt = Λ̂[1:end,(N_t+1):end]

Ω̂ = L̂_nt'*inv(L̂_nt*L̂_nt')*L̂_t

itt = Y_t3[3,1:end]' - Y_nt3[3,1:end]'*Ω̂

mean(itt)

quantile(itt', [0; 0.25; 0.5; 0.75; 1])

plot(density(itt'))

W_2t = Vector{Union{Missing, Float64}}(missing, 2165)

W_3t = Vector{Union{Missing, Float64}}(missing, 2165)

for i = 1:2165

for j = 6:10

if !ismissing(W_t[j,i])

if ismissing(W_2t[i])

W_2t[i] = W_t[j,i]

else

W_2t[i] = (W_2t[i] + W_t[j,i])/2

end

for j = 11:15

if !ismissing(W_t[j,i])

if ismissing(W_3t[i])

W_3t[i] = W_t[j,i]

else

W_3t[i] = (W_3t[i] + W_t[j,i])/2

end

inc = W_3t[index_t_na] - (Y_t3[3,1:end] - itt').*W_2t[index_t_na]./Y_t3[2,1:end]

quantile(skipmissing(inc), 0.5)

# 11.7