a final version for v 1.1

Project.toml

@@ -20,13 +20,13 @@ SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 [compat]
 julia = "1.0, 1.1, 1.2, 1.3"
 Combinatorics = "1"
-Distributions = "0.19, 0.20, 0.21, 0.22"
+Distributions = "0.19, 0.20, 0.21, 0.22, 0.23"
 QuadGK = "2"
-SpecialFunctions = "0.8, 0.9, 0.10"
-StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32"
+SpecialFunctions = "0.8, 0.9, 0.10, 0.11"
+StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33"
 HCubature = "1.4"
 HypothesisTests = "0.9"
-Roots = "0.7, 0.8"
+Roots = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/corgen.jl

@@ -44,12 +44,16 @@ end
   # generates covariance matrix
 
 """
-    cormatgen(n::Int, ρ::Float64 = 0.5, ordered::Bool = false, altersing::Bool = true)
+    cormatgen(n::Int = 20)
 
-Returns symmetric correlation matrix Σ of size n x n, with reference correlation 0 < ρ < 1.
-If ordered = false, Σ elements varies arround ρ, i.e. σᵢⱼ ≈ ρ+δ else they drop
-as indices differences rise, i.e. σᵢⱼ ≳ σᵢₖ as |i-j| < |i-k|.
-If altersing = true, some σ are positive and some negative, else ∀ᵢⱼ σᵢⱼ ≥ 0.
+Returns symmetric correlation matrix Σ of size n x n.
+Method:
+
+    a = rand(n,n)
+    b = a*a'
+    c = b./maximum(b)
+
+Example:
 
 ```jldoctest
 julia> Random.seed!(43);
@@ -68,6 +72,28 @@ function cormatgen(n::Int = 20)
 
 end
 
+"""
+    cormatgen_rand(n::Int = 20)
+
+Returns symmetric correlation matrix Σ of size n x n.
+Method:
+
+    a = rand(n,n)
+    b = a*a'
+    diagb = Matrix(Diagonal(1 ./sqrt.(LinearAlgebra.diag(b))))
+    b = diagb*b*diagb
+
+In general gives higher correlations than the cormatgen(). Example:
+
+```jldoctest
+julia> Random.seed!(43);
+
+julia> cormatgen_rand(2)
+2×2 Array{Float64,2}:
+ 1.0       0.879086
+ 0.879086  1.0
+```
+"""
 function cormatgen_rand(n::Int = 20)
   a = rand(n,n)
   b = a*a'
@@ -76,6 +102,18 @@ function cormatgen_rand(n::Int = 20)
   (b+b')/2
 end
 
+"""
+    cormatgen_constant(n::Int, α::Float64)
+
+Returns the constant correlation matrix with constant correlations equal to 0 <= α <= 1
+
+```julia
+julia> cormatgen_constant(2, 0.5)
+2×2 Array{Float64,2}:
+ 1.0  0.5
+ 0.5  1.0
+```
+"""
 function cormatgen_constant(n::Int, α::Float64)
   @assert 0 <= α <= 1 "α should satisfy 0 <= α <= 1"
   α .*ones(n, n) .+(1-α) .*Matrix(1.0I, n,n)
@@ -86,6 +124,22 @@ function random_unit_vector(dim::Int)
   result /= norm(result)
 end
 
+"""
+    cormatgen_constant_noised(n::Int, α::Float64; ϵ::Float64 = (1 .-α)/2.)
+
+Returns the constant correlation matrix of size n x n  with constant correlations equal to 0 <= α <= 1
+and additinal noise determinde by ϵ.
+
+```julia
+julia> Random.seed!(43);
+
+julia> cormatgen_constant_noised(3, 0.5)
+3×3 Array{Float64,2}:
+ 1.0       0.506271  0.285793
+ 0.506271  1.0       0.475609
+ 0.285793  0.475609  1.0
+```
+"""
 function cormatgen_constant_noised(n::Int, α::Float64; ϵ::Float64 = (1 .-α)/2.)
   @assert 0 <= ϵ <= 1-α "ϵ must satisfy 0 <= ϵ <= 1-α"
   result = cormatgen_constant(n, α)
@@ -94,6 +148,30 @@ function cormatgen_constant_noised(n::Int, α::Float64; ϵ::Float64 = (1 .-α)/2
   result - ϵ .*Matrix(1.0I, size(result)...)
 end
 
+"""
+    cormatgen_two_constant(n::Int, α::Float64, β::Float64)
+
+Returns the correlation matrix of size n x n of correlations determined by two constants, first must be greater than the second.
+
+```julia
+julia> cormatgen_two_constant(6, 0.5, 0.1)
+6×6 Array{Float64,2}:
+ 1.0  0.5  0.5  0.1  0.1  0.1
+ 0.5  1.0  0.5  0.1  0.1  0.1
+ 0.5  0.5  1.0  0.1  0.1  0.1
+ 0.1  0.1  0.1  1.0  0.1  0.1
+ 0.1  0.1  0.1  0.1  1.0  0.1
+ 0.1  0.1  0.1  0.1  0.1  1.0
+
+
+ julia> cormatgen_two_constant(4, 0.5, 0.1)
+ 4×4 Array{Float64,2}:
+  1.0  0.5  0.1  0.1
+  0.5  1.0  0.1  0.1
+  0.1  0.1  1.0  0.1
+  0.1  0.1  0.1  1.0
+```
+"""
 function cormatgen_two_constant(n::Int, α::Float64, β::Float64)
   @assert α > β "First argument must be greater"
   result = fill(β, (n,n))
@@ -102,6 +180,23 @@ function cormatgen_two_constant(n::Int, α::Float64, β::Float64)
   result
 end
 
+"""
+    cormatgen_two_constant_noised(n::Int, α::Float64, β::Float64; ϵ::Float64= (1-α)/2)
+
+Returns the correlation matrix of size n x n  of correlations determined by two constants, first must be greater than the second.
+Additional noise is introduced by the parameter ϵ.
+
+```julia
+julia> Random.seed!(43);
+
+julia> cormatgen_two_constant_noised(4, 0.5, 0.1)
+4×4 Array{Float64,2}:
+  1.0         0.314724   0.190368  -0.0530078
+  0.314724    1.0       -0.085744   0.112183
+  0.190368   -0.085744   1.0        0.138089
+ -0.0530078   0.112183   0.138089   1.0
+```
+"""
 function cormatgen_two_constant_noised(n::Int, α::Float64, β::Float64; ϵ::Float64= (1-α)/2)
   @assert ϵ <= 1-α
   result = cormatgen_two_constant(n, α, β)
@@ -110,11 +205,53 @@ function cormatgen_two_constant_noised(n::Int, α::Float64, β::Float64; ϵ::Flo
   result - ϵ .*Matrix(1.0I, size(result)...)
 end
 
+"""
+    cormatgen_toeplitz(n::Int, ρ::Float64)
+
+Returns the correlation matrix of size n x n of the Toeplitz structure with
+maximal correlation equal to ρ.
+
+```julia
+julia> cormatgen_toeplitz(5, 0.9)
+5×5 Array{Float64,2}:
+ 1.0     0.9    0.81  0.729  0.6561
+ 0.9     1.0    0.9   0.81   0.729
+ 0.81    0.9    1.0   0.9    0.81
+ 0.729   0.81   0.9   1.0    0.9
+ 0.6561  0.729  0.81  0.9    1.0
+
+julia> cormatgen_toeplitz(5, 0.6)
+5×5 Array{Float64,2}:
+ 1.0     0.6    0.36  0.216  0.1296
+ 0.6     1.0    0.6   0.36   0.216
+ 0.36    0.6    1.0   0.6    0.36
+ 0.216   0.36   0.6   1.0    0.6
+ 0.1296  0.216  0.36  0.6    1.0
+```
+"""
 function cormatgen_toeplitz(n::Int, ρ::Float64)
   @assert 0 <= ρ <= 1 "ρ needs to satisfy 0 <= ρ <= 1"
   [ρ^(abs(i-j)) for i=0:n-1, j=0:n-1]
 end
 
+"""
+    cormatgen_toeplitz_noised(n::Int, ρ::Float64; ϵ=(1-ρ)/(1+ρ)/2)
+
+Returns the correlation matrix of size n x n of the Toeplitz structure with
+maximal correlation equal to ρ. Additiona noise id added by the ϵ parameter.
+
+```julia
+julia> Random.seed!(43);
+
+julia> cormatgen_toeplitz_noised(5, 0.9)
+5×5 Array{Float64,2}:
+ 1.0       0.89656   0.812152  0.720547  0.64318
+ 0.89656   1.0       0.918136  0.832571  0.734564
+ 0.812152  0.918136  1.0       0.915888  0.822804
+ 0.720547  0.832571  0.915888  1.0       0.903819
+ 0.64318   0.734564  0.822804  0.903819  1.0
+```
+"""
 function cormatgen_toeplitz_noised(n::Int, ρ::Float64; ϵ=(1-ρ)/(1+ρ)/2)
   @assert 0 <= ϵ <= (1-ρ)/(1+ρ) "ϵ must satisfy 0 <= ϵ <= (1-ρ)/(1+ρ)"
 

test/test_on_data/corelations.jl

@@ -11,6 +11,7 @@ using PyCall
 using JLD2
 using FileIO
 using Distributed
+using HCubature
 
 #using QuadGK
 #using Cubature
@@ -26,7 +27,7 @@ addprocs(3)
 nprocs()
 
 
-function ccl(u::Vector, θ::Union{Int, Float64})
+function ccl(u, θ::Float64)
   q = length(u)
   (1-q+sum(u.^(-θ)))^(-q-1/θ)*prod([u[j]^(-θ-1)*((j-1)*θ+1) for j in 1:q])
 end
@@ -66,15 +67,14 @@ function c4symmarg(θ::Union{Int, Float64}, g::Function, d)
   [a,b,c,e,f]
 end
 
-function c3empirical(θ::Union{Float64, Int}, cop::String, d, rev = false, t::Int = 1_000_000)
-  u = quantile.(d, archcopulagen(t, 3, θ, cop; rev = rev))
+function c3empirical(cop, d, t::Int = 1_000_000)
+  u = quantile.(d, simulate_copula(t, cop))
   c = cumulants(u, 3)[3]
-  println(θ)
   [c[1,1,1], c[1,1,2], c[1,2,3]]
 end
 
-function c4empirical(θ::Union{Float64, Int}, cop::String, d, rev::Bool= false, t::Int = 5000000)
-  u = quantile.(d, archcopulagen(t, 4, θ, cop; rev = rev))
+function c4empirical(cop, d, t::Int = 5000000)
+  u = quantile.(d, simulate_copula(t, cop))
   c = cumulants(u, 4)[4]
   [c[1,1,1,1], c[1,1,1,2], c[1,1,2,2], c[1,1,2,3], c[1,2,3,4]]
 end
@@ -90,8 +90,8 @@ function stats3theoreunifcl()
 end
 
 
-function statstheoreticalclayton(m::Int = 3, d = Normal(0,1))
-  θ = map(i -> ρ2θ(i, "clayton"), 0.05:0.05:.95)
+function statstheoreticalclayton(m::Int = 3, d = Normal(0,1), ρ)
+  θ = map(i -> ρ2θ(i, "clayton"), ρ)
   n = length(θ)
   t = (m == 3) ? zeros(n, 3) : zeros(n, 5)
   for i in 1:n
@@ -103,31 +103,30 @@ end
 
 #save("c3clgmarg.jld2", "emp", e, "theoret", t, "rho", ρ)
 
-function empiricalcums(copula::String, m::Int, ρ::Vector{Float64}, rev = false)
-  θ = map(i -> ρ2θ(i, copula), ρ)
-  println(θ)
-  n = length(θ)
+function empiricalcums(copula, m::Int, ρ::Vector{Float64})
+  println(ρ)
+  n = length(ρ)
   e = (m == 3) ? zeros(n, 3) : zeros(n, 5)
   for i in 1:n
     println(i)
-    e[i,:] = (m == 3) ? c3empirical(θ[i], copula, Normal(0,1), rev) : c4empirical(θ[i], copula, Normal(0,1), rev)
+    e[i,:] = (m == 3) ? c3empirical(copula(3, ρ[i], "Spearman"), Normal(0,1)) : c4empirical(copula(4, ρ[i], "Spearman"), Normal(0,1))
   end
   e
 end
 
 
-function plotc3empth(e,t, ρ)
+function plotc3empth(e,t, ρ, ρ1)
   mpl.rc("font", family="serif", size = 7)
   fig, ax = subplots(figsize = (2.5, 2.))
   fx = matplotlib.ticker.ScalarFormatter()
   fx.set_powerlimits((-1, 4))
   ax.yaxis.set_major_formatter(fx)
   plot(ρ, e[:,3], "o", color = "black", label = "\$ \\mathbf{i} = (1,2,3) \$", markersize = 1.)
-  plot(ρ, t[:,3], color = "black", linewidth = 0.5)
+  plot(ρ1, t[:,3], color = "black", linewidth = 0.5)
   plot(ρ, e[:,2], "o", color = "blue", label = "\$ \\mathbf{i} = (1,1,2) \$", markersize = 1.)
-  plot(ρ, t[:,2], color = "blue", linewidth = 0.5)
+  plot(ρ1, t[:,2], color = "blue", linewidth = 0.5)
   plot(ρ, e[:,1], "o", color = "red", label = "\$ \\mathbf{i} = (1,1,1) \$", markersize = 1.)
-  plot(ρ, t[:,1], color = "red", linewidth = 0.5)
+  plot(ρ1, t[:,1], color = "red", linewidth = 0.5)
   PyPlot.ylabel("cumulant element", labelpad = -1.0)
   PyPlot.xlabel("Spearman \$ \\rho \$ between marginals", labelpad = -1.0)
   ax.legend(fontsize = 6, loc = 3, ncol = 1)
@@ -197,36 +196,48 @@ end
 
 if false
   ρ = collect(0.05:0.05:.95)
-  e4 = empiricalcums("gumbel", 4, ρ, false)
-  save("pics/c4gugmarg.jld2", "emp", e4,"rho", ρ)
-  t = statstheoreticalclayton(4, Normal(0,1));
+  e4 = empiricalcums(Clayton_cop, 4, ρ)
+  save("pics/c4clgmarg.jld2", "emp", e4,"rho", ρ)
+  t = statstheoreticalclayton(4, Normal(0,1), ρ);
   save("pics/teorc4clgmarg.jld2", "theor", t,"rho", ρ)
-
+end
+if false
+  ρ = collect(0.05:0.05:.95)
+  e3 = empiricalcums(Clayton_cop, 3, ρ)
+  save("pics/c3clgmarg.jld2", "emp", e3,"rho", ρ)
+  t = statstheoreticalclayton(3, Normal(0,1), ρ);
+  save("pics/teorc3clgmarg.jld2", "theor", t,"rho", ρ)
+end
+if false
   ρ = [i for i in 0.02:0.02:.85]
-  e = empiricalcums("frank", 3, ρ)
+  e = empiricalcums(Frank_cop, 3, ρ)
   save("pics/c3frgmarg.jld2", "emp", e,"rho", ρ)
-
+end
+if false
   ρ = [i for i in 0.02:0.02:.98]
-  e = empiricalcums("gumbel", 3, ρ)
+  e = empiricalcums(Gumbel_cop, 3, ρ)
   save("pics/c3gugmarg.jld2", "emp", e,"rho", ρ)
-
-  ρ = [i for i in 0.02:0.02:.5]
-  e = empiricalcums("amh", 3, ρ)
+end
+if false
+  ρ = vcat([i for i in 0.02:0.02:.49], [0.499])
+  e = empiricalcums(AMH_cop, 3, ρ)
   save("pics/c3amhgmarg.jld2", "emp", e,"rho", ρ)
 end
 
 
-c3els = load("pics/c3clgmarg.jld2")
-plotc3empth(c3els["emp"],c3els["theoret"], c3els["rho"])
+emp = load("pics/c3clgmarg.jld2")
+t = load("pics/teorc3clgmarg.jld2")
+plotc3empth(emp["emp"], t["theor"], emp["rho"], t["rho"])
 
 c3els = load("pics/c3gugmarg.jld2")
 plotc3emp(c3els["emp"], c3els["rho"], "gu")
 
+c3els = load("pics/c3frgmarg.jld2")
+plotc3emp(c3els["emp"], c3els["rho"], "fr")
+
 c3els = load("pics/c3amhgmarg.jld2")
 plotc3emp(c3els["emp"], c3els["rho"], "amh")
 
-c4els = load("pics/c4clgmarg.jld2")
-plotc4emp(c4els["emp"], c4els["rho"], "cl")
 
 emp = load("pics/c4clgmarg.jld2")
 t = load("pics/teorc4clgmarg.jld2")