Skip to content

Commit

Permalink
first step of new dispatching
Browse files Browse the repository at this point in the history
  • Loading branch information
@kdomino
kdomino committed Apr 17, 2020
1 parent ddd0e32 commit 65aa2cf
Show file tree
Hide file tree
Showing 5 changed files with 110 additions and 80 deletions.
1 change: 1 addition & 0 deletions src/DatagenCopulaBased.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module DatagenCopulaBased
export chainfrechetcopulagen
export gcop2tstudent, gcop2frechet, gcop2marshallolkin
export nested_gumbel, nested_clayton, nested_frank, nested_amh
export chain_frechet, chain_archimedeans, rev_chain_archimedeans

# obsolete implemntations
export tstudentcopulagen, gausscopulagen, frechetcopulagen, marshallolkincopulagen
Expand Down
64 changes: 42 additions & 22 deletions src/chaincopulagendat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,40 +35,38 @@ function rand2cop(u1::Vector{Float64}, θ::Union{Int, Float64}, copula::String)
end

"""
chaincopulagen(t::Int, θ::Union{Vector{Float64}, Vector{Int}}, copula::String;
rev::Bool = false, cor::String = "")
chain_archimedeans(t::Int, θ::Union{Vector{Float64}, Vector{Int}},
copula::Union{String, Vector{String}};
cor::String = "")
Returns: t x n Matrix{Float}, t realisations of n variate data, where n = length(θ)+1.
To generate data uses chain of bivariate Archimedean one parameter copulas.
Following copula families are supported: clayton, frank and amh -- Ali-Mikhail-Haq.
If rev == true, reverse the copula output i.e. u → 1-u (we call it reversed copula).
It cor == pearson, kendall, uses correlation coeficient as a parameter
```jldoctest
julia> Random.seed!(43);
julia> chaincopulagen(10, [4., 11.], "frank")
10×3 Array{Float64,2}:
0.180975 0.386303 0.879254
0.775377 0.247895 0.144803
0.888934 0.426854 0.772457
0.924876 0.395564 0.223155
0.408278 0.139002 0.142997
0.912603 0.901252 0.949828
0.828727 0.0295759 0.0897796
0.400537 0.0337673 0.27872
0.429437 0.462771 0.425435
0.955881 0.953623 0.969038
julia> chain_archimedeans(1, [4., 11.], "frank")
1×3 Array{Float64,2}:
0.180975 0.492923 0.679345
julia> Random.seed!(43);
julia> chain_archimedeans(1, [4., 11.], ["frank", "clayton"])
1×3 Array{Float64,2}:
0.180975 0.492923 0.600322
```
"""
VFI = Union{Vector{Float64}, Vector{Int}}
chaincopulagen(t::Int, θ::VFI, copula::String; rev::Bool = false, cor::String = "") =
chaincopulagen(t, θ, [fill(copula, length(θ))...] ; rev = rev, cor= cor)
chain_archimedeans(t::Int, θ::VFI, copula::String; cor::String = "") =
chain_archimedeans(t, θ, [fill(copula, length(θ))...]; cor = cor)


function chaincopulagen(t::Int, θ::VFI, copula::Vector{String}; rev::Bool = false, cor::String = "")
function chain_archimedeans(t::Int, θ::VFI, copula::Vector{String}; cor::String = "")
if ((cor == "Spearman") | (cor == "Kendall"))
θ = map(i -> usebivρ(θ[i], copula[i], cor), 1:length(θ))
else
Expand All @@ -78,7 +76,29 @@ function chaincopulagen(t::Int, θ::VFI, copula::Vector{String}; rev::Bool = fal
for i in 1:length(θ)
u = hcat(u, rand2cop(u[:, i], θ[i], copula[i]))
end
rev ? 1 .-u : u
return u
end

"""
rev_chain_archimedeans(t::Int, θ::VFI, copula::Union{String, Vector{String}};
cor::String = "")
Returns 1 .- chain_archimedeans(t, θ, copula; cor)
```jldoctest
julia> Random.seed!(43);
julia> rev_chain_archimedeans(1, [4., 11.], ["frank", "clayton"])
1×3 Array{Float64,2}:
0.819025 0.507077 0.399678
```
"""

function rev_chain_archimedeans(t::Int, θ::VFI, copula::Union{String, Vector{String}};
cor::String = "")
1 .- chain_archimedeans(t, θ, copula; cor = cor)
end

"""
Expand Down Expand Up @@ -123,15 +143,15 @@ end
# chain frechet copulas

"""
chainfrechetcopulagen(t::Int, α::Vector{Float64}, β::Vector{Float64} = zeros(α))
chain_frechet(t::Int, α::Vector{Float64}, β::Vector{Float64} = zeros(α))
Retenares data from a chain of bivariate two parameter frechet copuls with parameters
vectors α and β, such that ∀ᵢ 0 α[i] + β[i] ≤1 α[i] > 0, and β[i] > 0 |α| = |β|
```jldoctest
julia> Random.seed!(43)
julia> julia> chainfrechetcopulagen(10, [0.6, 0.4], [0.3, 0.5])
julia> chain_frechet(10, [0.6, 0.4], [0.3, 0.5])
10×3 Array{Float64,2}:
0.996764 0.996764 0.996764
0.204033 0.795967 0.204033
Expand All @@ -145,7 +165,7 @@ julia> julia> chainfrechetcopulagen(10, [0.6, 0.4], [0.3, 0.5])
0.804096 0.851275 0.955881
```
"""
function chainfrechetcopulagen(t::Int, α::Vector{Float64}, β::Vector{Float64} = zero(α))
function chain_frechet(t::Int, α::Vector{Float64}, β::Vector{Float64} = zero(α))
length(α) == length(β) || throw(AssertionError("different lengths of parameters"))
minimum(α) >= 0 || throw(DomainError("negative α parameter"))
minimum(β) >= 0 || throw(DomainError("negative β parameter"))
Expand Down
45 changes: 39 additions & 6 deletions src/copulagendat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,17 +6,24 @@
supports following copulas:
elliptical
- gaussian_cop,
args = Σ::Matrix{Float} (the correlation matrix, the covariance one will be normalised)
args: Σ::Matrix{Float} (the correlation matrix, the covariance one will be normalised)
- tstudent_cop,
args = Σ::Matrix{Float}, ν::Int
args: Σ::Matrix{Float}, ν::Int
- frechet,
args = n::Int, α::Union{Int, Float} (number of marginals, parameter of maximal copula α ∈ [0,1])
args: n::Int, α::Union{Int, Float} (number of marginals, parameter of maximal copula α ∈ [0,1])
or
args = n::Int, α::Union{Int, Float}, β::Union{Int, Float}, supported only for n = 2 and α, β, α+β ∈ [0,1]
args: n::Int, α::Union{Int, Float}, β::Union{Int, Float}, supported only for n = 2 and α, β, α+β ∈ [0,1]
- chain_frechet, simulates a chain of bivariate frechet copulas
args: α::Vector{Float64}, β::Vector{Float64} = zero(α) - vectors of
parameters of subsequent bivariate copulas. Each element must fulfill
conditions of the frechet copula. n.o. marginals = length(α)+1.
We require length(α) = length(β).
- marshalolkin,
args = λ::Vector{Float64} λ = [λ₁, λ₂, ..., λₙ, λ₁₂, λ₁₃, ..., λ₁ₙ, λ₂₃, ..., λₙ₋₁ₙ, λ₁₂₃, ..., λ₁₂...ₙ]
args: λ::Vector{Float64} λ = [λ₁, λ₂, ..., λₙ, λ₁₂, λ₁₃, ..., λ₁ₙ, λ₂₃, ..., λₙ₋₁ₙ, λ₁₂₃, ..., λ₁₂...ₙ]
n.o. margs = ceil(Int, log(2, length(λ)-1)), params:
- archimedean: gumbel, clayton, amh (Ali-Mikhail-Haq), frank.
Expand All @@ -34,8 +41,18 @@
- Ψ::Vector{Vector{Float64}}, ϕ::Vector{Float64}, θ::Float64 (params of ground childeren, children and parent copulas)
- hierarchical nested (only Gumbel)
- nested_gumbel
args: θ::Vector{Float64} - vector of parameters from ground grounf child to parent
args: θ::Vector{Float64} - vector of parameters from ground ground child to parent
all copuals are bivariate n = length(θ)+1
- chain_archimedeans simulate the chain of bivariate archimedean copula,
args: - θ::Union{Vector{Float64}, Vector{Int}} - parameters of subsequent bivariate copulas
- copula::Union{Vector{String}, String}, indicates a bivariate copulas
or their sequence, supported are supports: clayton, frank and amh famillies
- keyword cor, if cor = ["Spearman", "Kendall"] uses these correlations in place of parameters
of subsequent buvariate copulas
- rev_chain_archimedeans reversed version of the chain_archimedeans
"""

function simulate_copula(t::Int, copula::Function, args...; cor = "")
Expand Down Expand Up @@ -114,3 +131,19 @@ function nestedarchcopulagen(t::Int, θ::Vector{Float64}, copula::String = "gumb
copula == "gumbel" || throw(AssertionError("generator supported only for gumbel familly"))
simulate_copula(t, nested_gumbel, θ)
end


function chainfrechetcopulagen(t::Int, α::Vector{Float64}, β::Vector{Float64} = zero(α))
simulate_copula(t, chain_frechet, α, β)
end

VFI = Union{Vector{Float64}, Vector{Int}}

function chaincopulagen(t::Int, θ::VFI, copula::Union{Vector{String}, String};
rev::Bool = false, cor::String = "")
if rev == false
return simulate_copula(t, chain_archimedeans, θ, copula, cor = cor)
else
return simulate_copula(t, rev_chain_archimedeans, θ, copula, cor = cor)
end
end
46 changes: 0 additions & 46 deletions src/nestedarchcopulagendat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,49 +228,3 @@ function testnestedθϕ(n::Vector{Int}, ϕ::Vector{Float64}, θ::Float64, copula
length(n) == length(ϕ) || throw(AssertionError("number of subcopulas ≠ number of parameters"))
(copula != "clayton") | (maximum(ϕ) < θ+2*θ^2+750*θ^5) || warn("θ << ϕ for clayton nested copula, marginals may not be uniform")
end


#=
function nestedarchcopulagen1(t::Int, n::Vector{Vector{Int}}, Ψ::Vector{Vector{Float64}},
ϕ::Vector{Float64}, θ::Float64,
copula::String = "gumbel")
copula == "gumbel" ||
throw(AssertionError("double nested cop. generator supported only for gumbel familly"))
θ <= minimum(ϕ) || throw(DomainError("wrong heirarchy of parameters"))
X = nested_gumbel(t, n[1], Ψ[1], ϕ[1]./θ)
for i in 2:length(n)
X = hcat(X, nested_gumbel(t, n[i], Ψ[i], ϕ[i]./θ))
end
phi(-log.(X)./levygen(θ, rand(t)), θ, "gumbel")
end
function nestedarchcopulagen1(t::Int, θ::Vector{Float64}, copula::String = "gumbel")
copula == "gumbel" ||
throw(AssertionError("hierarchically nasted cop. generator supported only for gumbel familly"))
testθ(θ[end], "gumbel")
issorted(θ; rev=true) || throw(DomainError("wrong heirarchy of parameters"))
θ = vcat(θ, [1.])
X = rand(t,1)
for i in 1:length(θ)-1
X = nestedstep("gumbel", hcat(X, rand(t)), ones(t), θ[i], θ[i+1])
end
X
end
function nestedarchcopulagen1(t::Int, n::Vector{Int}, ϕ::Vector{Float64}, θ::Float64, copula::String, m::Int = 0)
n1 = vcat([collect(1:n[1])], [collect(cumsum(n)[i]+1:cumsum(n)[i+1]) for i in 1:length(n)-1])
nestedarchcopulagen(t, sum(n)+m, n1, ϕ, θ, copula)
end
function nestedarchcopulagen1(t::Int, n::Int, n1::Vector{Vector{Int}}, ϕ::Vector{Float64}, θ::Float64,
copula::String)
copula in ["clayton", "amh", "frank", "gumbel"] || throw(AssertionError("$(copula) copula is not supported"))
nestedcopulag(copula, n1, ϕ, θ, rand(t, n+1))
end
=#
34 changes: 28 additions & 6 deletions test/chaincopulastests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,20 @@ end

@testset "chain of Archimedean copulas" begin
Random.seed!(43)
@test chain_archimedeans(1, [4., 11.], "frank") ≈ [0.180975 0.492923 0.679345] atol=1.0e-5
Random.seed!(43)
@test chain_archimedeans(1, [4., 11.], ["frank", "clayton"]) ≈ [0.180975 0.492923 0.600322] atol=1.0e-5
Random.seed!(43)
@test rev_chain_archimedeans(1, [4., 11.], ["frank", "clayton"]) ≈ [0.819025 0.507077 0.399678] atol=1.0e-5

cops = ["clayton", "clayton", "clayton", "frank", "amh", "amh"]
x = chaincopulagen(500000, [-0.9, 3., 2, 4., -0.3, 1.], cops)
Random.seed!(43)
x = simulate_copula(500000, chain_archimedeans, [-0.9, 3., 2, 4., -0.3, 1.], cops)
#test obsolete implemntation
Random.seed!(43)
x1 = chaincopulagen(500000, [-0.9, 3., 2, 4., -0.3, 1.], cops)
@test norm(x-x1) ≈ 0.

@test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
@test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
@test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
Expand All @@ -39,17 +51,21 @@ end
end
@testset "correlations" begin
Random.seed!(43)
x = chaincopulagen(500000, [0.6, -0.2], "clayton"; cor = "Spearman")
x = simulate_copula(500000, chain_archimedeans, [0.6, -0.2], "clayton"; cor = "Spearman")
@test corspearman(x[:,1], x[:,2]) ≈ 0.6 atol=1.0e-2
@test corspearman(x[:,2], x[:,3]) ≈ -0.2 atol=1.0e-2
Random.seed!(43)
x = chaincopulagen(500000, [0.6, -0.2], "clayton"; cor = "Kendall")
x = simulate_copula(500000, chain_archimedeans, [0.6, -0.2], "clayton"; cor = "Kendall")
@test corkendall(x[:,1], x[:,2]) ≈ 0.6 atol=1.0e-3
@test corkendall(x[:,2], x[:,3]) ≈ -0.2 atol=1.0e-3
end
@testset "rev copula" begin
Random.seed!(43)
x = chaincopulagen(500000, [-0.9, 3., 2.], "clayton"; rev = true)
x = simulate_copula(500000, rev_chain_archimedeans, [-0.9, 3., 2.], "clayton")
# test obsolete dispatching
Random.seed!(43)
x1 = chaincopulagen(500000, [-0.9, 3., 2.], "clayton"; rev = true)
@test norm(x - x1) ≈ 0.
@test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
@test tail(x[:,3], x[:,4], "r") ≈ 1/(2^(1/2)) atol=1.0e-1
@test tail(x[:,3], x[:,4], "l", 0.0001) ≈ 0
Expand All @@ -58,7 +74,13 @@ end
@testset "chain of Frechet copulas" begin
@test fncopulagen([0.2, 0.4], [0.1, 0.1], [0.2 0.4 0.6; 0.3 0.5 0.7]) == [0.6 0.4 0.2; 0.7 0.5 0.3]
Random.seed!(43)
x = chainfrechetcopulagen(500000, [0.9, 0.6, 0.2])
@test chain_frechet(1, [0.6, 0.4], [0.3, 0.5]) ≈ [0.888934 0.775377 0.180975] atol=1.0e-5
Random.seed!(43)
x = simulate_copula(500000, chain_frechet, [0.9, 0.6, 0.2])
# test obsolete dispatching
Random.seed!(43)
x1 = chainfrechetcopulagen(500000, [0.9, 0.6, 0.2])
@test norm(x - x1) ≈ 0.
@test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
@test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
@test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
Expand All @@ -67,7 +89,7 @@ end
@test tail(x[:,1], x[:,2], "l") ≈ 0.9 atol=1.0e-1
@test tail(x[:,1], x[:,4], "r") ≈ 0.2 atol=1.0e-1
Random.seed!(43)
x = chainfrechetcopulagen(500000, [0.8, 0.5], [0.2, 0.3]);
x = simulate_copula(500000, chain_frechet, [0.8, 0.5], [0.2, 0.3]);
@test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
@test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
@test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
Expand Down

0 comments on commit 65aa2cf

Please sign in to comment.