Extreme value datasets for GS 240

extreme_data.jl

srand(2000) # make sure we get the same data

#############
# LOGNORMAL #
#############

using SpecialFunctions

# CDF and inverse CDF of standard normal
ϕ(x) = erfc(-x/√2) / 2
ϕinv(p) = -√2*erfcinv(2p)

# generate N samples from log-normal
function lognormal(μ, σ², N)
    p = rand(N)              # draw uniformly in [0,1)
    logx = μ + ϕinv.(p)*√σ²
    exp.(logx)
end

data₁ = lognormal(2, 2, 500)

#################
# LOGHYPERBOLIC
#################

using FastGaussQuadrature
using Interpolations

# PDF of log-hyperbolic
function f(x; α=1., ϕ=1., μ=1., δ=1.)
    C = 1. / ((1./ϕ + 1./α) * δ * √(ϕ*α) * x * besselk(1., δ * √(ϕ*α)))
    A = (ϕ+α)/2
    B = (ϕ-α)/2
    z = (log.(x) - μ)
    C * exp.(-A*√(δ^2 + z.^2) + B*z)
end

# case of interest
h(x) = f(x, α=1./.9, ϕ=10., μ=.3, δ=.6)

# quadrature points and weights
xs, ws = gausslegendre(10000)

# integrate f from a to b
function integrate(f, a, b)
    c, d = (b-a)/2, (b+a)/2
    c * ws ⋅ f.(c*xs + d)
end

# corresponding CDF
H(x) = integrate(h, 0., x)

# invert function F in [a,b] using N interpolation points
function invert(F, a, b, N)
    knots = linspace(a, b, N)
    values = [F(x) for x in knots]
    idx = sortperm(values)
    knots = knots[idx]
    values = values[idx]
    itp = interpolate((values,), knots, Gridded(Linear()))
    p -> itp[p]
end

# quantile function
Q = invert(H, 1e-4, 1e3, 2000)

data₂ = [Q(p) for p in rand(500)]

##########
# PARETO
##########

using Distributions

data₃ = rand(Pareto(.9, 10.), 150)

##########################################################

# save data to disk
for (i,data) in enumerate([data₁, data₂, data₃])
    writedlm("data$i.dat", data)
end