ECO539B - Homework 1

hw01.jl

using Random, Distributions
using ProgressMeter
using LaTeXTabulars, LaTeXStrings
Random.seed!(0303)


"""
Generate `M` simulation draws for `k` parameters for a design given by `θfactor`
Inputs:
    - θfactor ... scalar with factor for θ_i (the "design")
    - k ... scalar for number of parameters / 1st dimension of θ
    - M ... scalar for total number of simulations / 2nd dimension of θ
Output:
    - θ ... kxM array with simulated draws of θ_i
"""
function genθ(θfactor; k = 10, M = 5000)
    θ = Array{Float64}(undef, k, M)
    for i in 1:k
        d = Normal(i*θfactor, 1)
        θ[i, :] = rand(d, M)
    end
    return θ
end

θ1 = genθ(1)
θ2 = genθ(10)
θ3 = genθ(0.1)


"""
Compute the rank of a `k`-dimensional vector of `θ` parameters
Inputs:
    - θ_m ... kx1 array with simulated draws of `θ_i`
Output:
    - R_m ... kx1 rank vector for simulation draw `m`
"""
function computeRankVector(θ_m)
    k = size(θ_m)[1]
    R_m = zeros(Float64, k)
    for i in 1:k
        R_m[i] = 1.0 + sum(θ_m .< θ_m[i]) + 0.5 * (sum(θ_m .== θ_m[i]) - 1.0)
    end
    return R_m
end

# R1 = mapslices(computeRankVector, θ1; dims = 1)
# R2 = mapslices(computeRankVector, θ2; dims = 1)
# R3 = mapslices(computeRankVector, θ3; dims = 1)

function bootstrapRankArray(θ_m; N = 1000, centered = true)
    k = size(θ_m)[1]
    θestar = θ_m .+ randn(k, N)    # k x N
    R_bs = Array{Float64}(undef, k, N)
    for b in 1:N
        if centered
            R_bs[:, b] = computeRankVector(θestar[:, b]) - computeRankVector(θ_m)
        else
            R_bs[:, b] = computeRankVector(θestar[:, b])
        end
    end
    return sort(R_bs, dims = 2)
end

function bootstrapIntervals(θ_m; α = 0.05, N = 1000)
    Rjat = computeRankVector(θ_m)
    R = bootstrapRankArray(θ_m)
    lb = Int(N*α/2)
    ub = Int(N*(1-α/2))
    CI_ptile = [ Rjat - R[:, ub] Rjat - R[:, lb] ]
    CI_efron = [ Rjat + R[:, lb] Rjat - R[:, ub] ]
    return CI_ptile, CI_efron
end


function computeCoverageVector(ci)
    k = size(ci)[1]
    coverage = Array{Bool}(undef, size(ci)[1])
    for i in 1:k
        coverage[i] = (ci[i, 1] <= i && ci[i, 2] >= i)
    end
    return coverage
end


"""
Simulate coverage over M simulations for the percentile and Efron method
Inputs:
    - θ ... kxM array with simulated draws of θ_i
Output:
    - coverage_ptile, coverage_efron ... kx1 vectors with simulated coverage
"""
function simulateCoverage(θ)
    k, M  = size(θ)
    coverageArray_ptile = Array{Bool}(undef, k, M)
    coverageArray_efron = Array{Bool}(undef, k, M)
    @showprogress for m in 1:M
        θ_m = θ[:, m]
        ci_ptile, ci_efron = bootstrapIntervals(θ_m)
        coverageArray_ptile[:, m] = computeCoverageVector(ci_ptile)
        coverageArray_efron[:, m] = computeCoverageVector(ci_efron)
    end
    coverage_ptile = mean(coverageArray_ptile, dims = 2)
    coverage_efron = mean(coverageArray_efron, dims = 2)
    return coverage_ptile, coverage_efron
end

efron1, ptile1 = simulateCoverage(θ1)
efron2, ptile2 = simulateCoverage(θ2)
efron3, ptile3 = simulateCoverage(θ3)

results = [1:10 ptile1 ptile2 ptile3 efron1 efron2 efron3]

latex_tabular(
    "01-Bootstrap/hw01/coverage.tex",
    Tabular("rcccccc"),
    [Rule(:top),
    ["" , MultiColumn(3, :c, "Percentile"), MultiColumn(3, :c, "Efron")],
    CMidRule("lr", 2, 4), CMidRule("lr", 5, 7),
    [L"i", L"\theta_i = i", L"\theta_i = 10i", L"\theta_i = 0.1i", L"\theta_i = i", L"\theta_i = 10i", L"\theta_i = 0.1i"],
    Rule(:mid),
    results,
    Rule(:bottom)])