VerticalBargaining.jl

#%% Load libraries
using CSV, DataFrames, DataFramesMeta
using LinearAlgebra
using Optim # GMM() optimization
# using FixedEffectModels
# using RegressionTables
using StatsBase # sample() for bootstrap sample
using ProgressMeter
using Latexify

#%% Read data
cd("/home/alvaro/Dropbox/Princeton/2020-Spring/542-IO/PS2")
df = CSV.read("data/PS2_data.csv", header = [:s; :p; :x; :z]);

#%% Exogenous parameters
α  = -0.5
β₀ =  1.5
β₁ =  0.5

# θ = [ϕ; λ; γ₀; γ₁]

#%% Compute variables that come from the model
# u  = log.(df.s) .- log.(1 .- df.s);
df = @transform(df, ξ = - β₀ .- df.x .* β₁ .- α .* df.p .- log.( (1 .- df.s) ./ df.s))
df = @transform(df, ones = ones(size(df)[1]))
#df = @transform(df, sinv = 1 .+ 1 ./ α .* (1 .- :s));
# df = @transform(df, k = 1 ./ α .* (1 .- :s))
# df = @transform(df, p₂ = df.p .+ df.k)
#ω(λ, γ₁, ϕ) = α .* (1 .- df.s) .* (2 .- df.s) / ((1 .- df.s) .* α) .+ df.p .- λ .- γ₁ .* df.z .- ϕ


#%% Define GMM functions
function GMM(θ, α, df::DataFrame, Z::Array{Symbol,1})
    lengthZ = length(Z)
    g = zeros(lengthZ)
    N = size(df)[1]
    λ_γ = θ[1]
    ϕ   = θ[2]
    γ₁  = θ[3]
    # Loop through Z
    for i in 1:lengthZ
        z = Z[i]
        g[i] = (
            1 / N .* vec(df[:, z])' *
            ((df.p .- λ_γ .- df.z .* γ₁ .+ 1 ./ (α .* (1 .- df.s)))
            ./
            ((df.p .- λ_γ .- df.z .* γ₁ .+ 1 ./ (α .* (1 .- df.s))) .* df.s .- 1 ./  (α .* (1 .- df.s)))
            .- ϕ
            )
        )[1]
    end
    # Return object
    obj = 1 / N .* g' * g
    return obj[1]
end

multieqGMM(θ, α, df, Z)

function multieqGMM(initial, α, df::DataFrame, Z::Array{Symbol,1})
    Optim.minimizer(Optim.optimize(θ -> GMM(θ, α, df, Z), initial[1:3], LBFGS()))
end

#%% Run GMM to get point estimates
# Define initial parameter values:
ϕ = 0.5
λ = 1
λ_γ = 1
γ₁ = 0
θ  = [λ_γ; ϕ; γ₁]



#%% Define GMM function that computes bootstrap standard errors
function bootstrapGMM(K, α, df::DataFrame, Z::Array{Symbol,1}, initial=zeros(size(Z)[1]))
    coefs_data = multieqGMM(initial, α, df, Z)
    coefs_bootstrap = zeros(K, 3)
    @showprogress 1 for i in 1:K
        sample                = df[StatsBase.sample(axes(df, 1), size(df)[1];
                                   replace = true, ordered = false), :]
        coefs_bootstrap[i, :] = multieqGMM(initial, α, sample, Z)
    end
    GMMestimates = DataFrame()
    GMMestimates.params = ["λ+γ₀"; "ϕ"; "γ₁"]
    GMMestimates.coefs = vec(coefs_data)
    GMMestimates.SE = vec(StatsBase.std(coefs_bootstrap, dims = 1))
    return GMMestimates
end


#%% Compute estimates
# Number of bootstrap repetitions:
K = 10000

#
Z = [:ones, :x, :z, :ξ];
est1 = bootstrapGMM(K, α, df, Z)

Z = [:ones, :x, :z];
est2 = bootstrapGMM(K, α, df, Z)

Z = [:ones, :x, :ξ];
est3 = bootstrapGMM(K, α, df, Z)




#%% Create combined array of results to export LaTeX table
results = string.(zeros(6, 3))
for (index, est) in enumerate([est1, est2, est3])
    iter = 0
    for i in 1:3
        for j in 2:3
            iter = iter + 1
            if j == 2
                results[iter, index] = string.(round.(est[i, j], sigdigits = 3))
            else
                results[iter, index] = string.("(", round.(est[i, j], sigdigits = 3), ")")
            end
        end
    end
end
results

# Round results and convert to string for exporting
results = hcat(["λ+γ₀"; ""; "ϕ"; ""; "γ₁"; ""], results)
latexify(results, env = :tabular, head = [""; "(1)"; "(2)"; "(3)"], fmt = FancyNumberFormatter(3))