chapter5.jl

# Chapter 5

#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

# 5.1

# 5.2

# 5.2.1

# 5.2.2

rng = MersenneTwister(123457689)

N = 1000

υ = randn!(rng, zeros(N)).*3 .+ 1

p = 2

mean(υ .- p .> 0)

1 - cdf(Normal(1, 3), p)

# 5.2.3

p = rand!(rng, zeros(9)).*20 .- 10

s = Vector{Union{Missing, Float64}}(missing, 9)

for i = 1:9

s[i] = mean(υ .- p[i] .> 0)

end

a = ecdf(υ)

scatter(p, 1 .- s)

scatter!(ecdf(υ))

# note the R code has survival function rather than the cdf.

# 5.3.1

# 5.3.2

N = 1000

a = fill(2, N)

b = -3

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

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

y_star = a + b.*x + υ

y = y_star .> 0

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

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

scatter(x, y)

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

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

# 5.3.3

# 5.4

# 5.4.1

N = 100

p = 0.367

Head = rand!(rng, zeros(N)) .< p

mean(Head)

factorial(5)/(factorial(3)*factorial(2))

factorial(big(100))/(factorial(big(34))*factorial(big(66)))*(0.5.^100)

# 5.4.2

binom = function(p, N, p̂)

return (p̂*N*log(p) + (1 - p̂)*N*log(1 - p))

end

binom_wrap = function(p)

return -binom(p, N, mean(Head))

end

a = optimize(binom_wrap, 0, 1)

a.minimum

a.minimizer

p = 1:1000

log_lik = binom.(p./1000, N, mean(Head))

scatter(p./1000, log_lik)

vline!([mean(Head), 0.367])

# 5.4.3

# 5.4.4

# note that the codes is quite a bit different from the R code.

# standard normal density

ϕ = function(z)

(exp.(-z.^2))./√π

end

# log standard normal density

logϕ = function(z)

-z.^2 - √π

end

ols_ml = function(par, y, X)

σ = exp(par[1])

β = par[2:3]

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

return -sum(logϕ.(z) .- log(σ))

end

ols_ml_wrap = function(par)

return ols_ml(par, y, X)

end

global y = y_star

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

init = [0.0; 2.0; -3.0]

a = optimize(ols_ml_wrap, init)

exp(a.minimizer[1])

a.minimizer[2:3]

# 5.4.5

# 5.4.6

Φ = function(v)

cdf(Normal(), v)

end

logΦ = function(v)

logcdf(Normal(), v)

end

probit = function(β, y, X)

return (1 .- y).*logΦ.(-X*β) + y.*logΦ.(X*β)

end

probit_wrap = function(β)

return -sum(probit(β, y, X))

end

global y = data1.y

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

a = optimize(probit_wrap, [2.0, -3.0])

a.minimizer

# 5.4.7

glm1 = glm(@formula(y ~ x), data1, Binomial(), ProbitLink())

# 5.5.

# 5.5.1

# 5.5.2

# 5.5.3

Random.seed!(123456789)

N = 5000

X_A = hcat(ones(N), rand(rng, N), rand(rng, N))

X_B = hcat(ones(N), rand(rng, N), rand(rng, N))

β = [1; -2; 3]

# Probit

υ_A = randn(rng, N)

y = (X_A .- X_B)*β + υ_A .> 0

data1 = DataFrame(y = y, x₁=X_A[1:end,2] .- X_B[1:end,2], x₂=X_A[1:end,3] .- X_B[1:end,3])

glm1 = glm(@formula(y ~ x₁ + x₂), data1, Binomial(), ProbitLink())

# Logit

υ_A = rand.(Gumbel(), N)

υ_B = rand.(Gumbel(), N)

y = (X_A .- X_B)*β .+ υ_A .- υ_B .> 0

data2 = DataFrame(y = y, x₁=X_A[1:end,2] .- X_B[1:end,2], x₂=X_A[1:end,3] .- X_B[1:end,3])

glm2 = glm(@formula(y ~ x₁ + x₂), data2, Binomial(), LogitLink())

# 5.6

# 5.6.1

# 5.6.2

# 5.6.3

biprobit = function(β, ρ, y, W₁, W₂, K)

rng = MersenneTwister(123457689)

ϵ = 1e-10

N = size(y)[1]

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

p_00 = zeros(N)

p_10 = zeros(N)

υ_k0 = zeros(N)

for k = 1:K

υ_k0 = υ[k] .< -W₁*β

p_00 = p_00 .+ υ_k0.*Φ.((-W₂*β .- ρ*υ[k])./√(1 - ρ^2))

p_10 = p_10 .+ (1 .- υ_k0).*Φ.(((W₁ - W₂)*β .+ (1 - ρ)*υ[k])./√(1 - ρ^2))

end

loglik = (y[1:end,1] .== 0).*(y[1:end,2] .== 0).*log.(p_00./K .+ ϵ) + (y[1:end,1] .== 1).*(y[1:end,2] .== 0).*log.(p_10./K .+ ϵ) + (y[1:end,1] .== 0).*(y[1:end,2] .== 1).*log.(1 .- p_00./K .- p_10./K .+ ϵ)

end

biprobit_wrap = function(par)

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

β = par[2:end]

-sum(biprobit(β, ρ, y, W₁, W₂, K))

end

# 5.6.4

# 5.6.5

f_logit = function(β, y, W₁, W₂)

ϵ = 1e-20

p_10 = exp.(W₁*β)./(1 .+ exp.(W₁*β) .+ exp.(W₂*β))

p_01 = exp.(W₂*β)./(1 .+ exp.(W₁*β) .+ exp.(W₂*β))

p_00 = 1 .- p_10 .- p_01

loglik = (y[1:end,1] .== 0).*(y[1:end,2] .== 0).*log.(p_00 .+ ϵ) + (y[1:end,1] .== 1).*(y[1:end,2] .== 0).*log.(p_10 .+ ϵ) + (y[1:end,1] .== 0).*(y[1:end,2] .== 1).*log.(p_01 .+ ϵ)

end

f_logit_wrap = function(β)

-sum(f_logit(β, y, W₁, W₂))

end

# 5.6.6

Random.seed!(123456789)

N = 1000

μ = [0; 0]

ρ = 0.1

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

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

X₁ = hcat(rand!(rng, zeros(N)), rand!(rng, zeros(N)))

X₂ = hcat(rand!(rng, zeros(N)), rand!(rng, zeros(N)))

a = fill(-1, N, 2)

b = -3

c = 4

U = a + b.*X₁ + c.*X₂ + Υ

global y = hcat((U[1:end,1] .> 0).*(U[1:end,1] .> U[1:end,2]), (U[1:end,2] .> 0).*(U[1:end,2] .> U[1:end,1]))

init = [log(0.1); -1; -3; 4]

global W₁ = hcat(ones(N), X₁[1:end,1], X₂[1:end,1])

global W₂ = hcat(ones(N), X₁[1:end,2], X₂[1:end,2])

global K = 100

a₁ = optimize(biprobit_wrap, init)

a₁.minimizer[2:end]

exp(a₁.minimizer[1])/(1 + exp(a₁.minimizer[1]))

b₁ = optimize(f_logit_wrap, init[2:end])

b₁.minimizer

# 5.7

# 5.7.1

path = "ENTER PATH HERE"

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

choice = Vector{Union{Missing, String}}(missing, 129696)

home = Vector{Union{Missing, Float64}}(missing, 129696)

income = Vector{Union{Missing, Float64}}(missing, 129696)

density = Vector{Union{Missing, Float64}}(missing, 129696)

index = []

for i = 1:129696

if x.CAR[i] == 1

choice[i] = "car"

end

if x.BUS[i] == 1

choice[i] = "bus"

end

if x.TRAIN[i] == 1

choice[i] = "train"

end

if x.HOMEOWN[i] == 1

home[i] = 1

elseif x.HOMEOWN[i] > 1

home[i] = 0

end

if x.HHFAMINC[i] > 0

income[i] = x.HHFAMINC[i]

end

if x.HTPPOPDN[i] > 0

density[i] = x.HTPPOPDN[i]/1000

end

if x.WRKCOUNT[i] > 0 && (x.MSACAT[i] == 1 || x.MSACAT[i] == 2)

index = [index; i]

end

#println(i)

end

car1 = choice .== "car"

train1 = choice .== "train"

urban1 = x.URBAN .== 1

rail = x.RAIL .== 1

x1 = DataFrame(car1=car1[index], train1=train1[index], home=home[index], HHSIZE=x.HHSIZE[index], income=income[index], urban1=urban1[index], density=density[index], MSACAT=x.MSACAT[index], rail=rail[index])

x1 = dropmissing(x1)

# Note in the R code this is y, which leads to confusion below.

vars = ["car1", "train1", "home", "HHSIZE", "income", "urban1", "density"]

summ_tab = Array{Union{Missing,Float64,String}}(missing,length(vars)+1,3)

summ_tab[1,2:end] = ["Rail" "No Rail"]

summ_tab[2:end,1] = vars

for i in 1:length(vars)

summ_tab[i+1,2] = mean(x1[1:end,i][x1.rail .== 1])

summ_tab[i+1,3] = mean(x1[1:end,i][x1.rail .== 0])

end

summ_tab

# 5.7.2

# without rail

x1_nr = x1[x1.rail .== 0,1:end]

glm_nr = glm(@formula(car1 ~ home + HHSIZE + income + urban1 + density), x1_nr, Binomial())

# with rail

x1_r = x1[x1.rail .== 1,1:end]

glm_r = glm(@formula(car1 ~ home + HHSIZE + income + urban1 + density), x1_r, Binomial())

# 5.7.3

Nᵣ = size(x1_r)[1]

global W₁ = hcat(ones(Nᵣ), x1_r.home, x1_r.HHSIZE, x1_r.income, x1_r.urban1, zeros(Nᵣ))

global W₂ = hcat(ones(Nᵣ), zeros(Nᵣ), zeros(Nᵣ), zeros(Nᵣ), zeros(Nᵣ), x1_r.density)

global y = hcat(x1_r.car1,x1_r.train1)

global K = 200

init = coef(glm_r)

a₂ = optimize(f_logit_wrap, init)

a₂.minimizer

init = [0; a₂.minimizer]

a₁ = optimize(biprobit_wrap, init, show_trace=true, iterations=5000)

a₁.minimizer[2:end]

exp(a₁.minimizer[1])/(1 + exp(a₁.minimizer[1]))

# 5.7.4

N_nr = size(x1_nr)[1]

W₁ = hcat(ones(N_nr), x1_nr.home, x1_nr.HHSIZE, x1_nr.income, x1_nr.urban1, zeros(N_nr))

W₂ = hcat(ones(N_nr), zeros(N_nr), zeros(N_nr), zeros(N_nr), zeros(N_nr), x1_nr.density)

# Probit Estimate

rng = MersenneTwister(123457689)

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

β̂ = a₁.minimizer[2:end]

K = 200

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

p_00 = zeros(N_nr)

p_10 = zeros(N_nr)

for k = 1:K

υ_k0 = υ[k] .< -W₁*β̂

p_00 = p_00 .+ υ_k0.*Φ.((-W₂*β̂ .- ρ̂*υ[k])./√(1 - ρ̂^2))

p_10 = p_10 .+ (1 .- υ_k0).*Φ.(((W₁ - W₂)*β̂ .+ (1 - ρ̂)*υ[k])./√(1 - ρ̂^2))

end

p_train_pr = (K .- p_00 .- p_10)./K

mean(p_train_pr)

# Logit Estimate

β̂ = a₂.minimizer

p_train_lg = exp.(W₂*β̂)./(1 .+ exp.(W₁*β̂) + exp.(W₂*β̂))

mean(p_train_lg)

# 5.8