Embedded Files
# 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)).*σ
Ȳ = F*Λ
Ys = Ȳ + 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
ŷ_ols = hcat(ones(5), Ys[31:35,2:7])*ω̂
res_ols = √mean((Ys[31:35,1] - ŷ_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))
ŷ_sc = Ys[31:35,2:end]*ω̂
res_sc = √mean((Ys[31:35,1] - ŷ_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
end
return ω_th
end
ω̂ = lasso(y_r, y_t, Y_r, Y_t, K=5)
X = hcat(ones(5), Ys[31:35,2:end])
ŷ_lasso = X*ω̂
res_lasso = √mean((Ys[31:35,1] - ŷ_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]
ŷ_fm = Yⱼ*Λ̂ⱼ'*inv(Λ̂ⱼ*Λ̂ⱼ')*λ̂ᵢ
res_fm = √mean((Ys[31:35,1] - ŷ_fm).^2)
plot(Ys[31:35,1],Ys[31:35,1], label = :none)
scatter!(Ys[31:35,1], ŷ_ols, label = "OLS")
scatter!(Ys[31:35,1], ŷ_sc, label = "SC")
scatter!(Ys[31:35,1], ŷ_lasso, label = "LASSO")
scatter!(Ys[31:35,1], ŷ_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
end
end
if x_names[1][j] == inc_name
for k = 1:8984
if x[k,j] >= 0
W[i,k] = x[k,j]
end
end
end
end
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
end
if !ismissing(W[11,i])
if W[11,i] == 0
treat[i] = 1
end
end
if !ismissing(Y[11,i])
if Y[11,i] == 0
treat[i] = 1
end
end
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
end
if length(X) > 1
treat = X[1:end,1]
for i = 2:N
treat = vcat(treat, X[1:end,i])
end
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
end
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
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
end
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
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
end
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
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
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
end
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
end
end
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
Google Sites
Report abuse