copem2D.jl

## Module `copem2D` for calculations related to bivariate empirical copulas
## Author: Arturo Erdely
## Version date: 2024-03-30
"""
# Module: `copem2D`
Given an observed random sample from a continuous bivariate random vector (X,Y)
this module provides functions to calculate:

- summary marginal and dependence statistics
- marginal and bivariate plots

> Dependence on external packages: `StatsPlots`

> Dependence on external files: `copem2Daux.jl`

Main functions: `dstat` `dplot` `resumen`
"""
module copem2D

export resumen, dstat, dplot 

using Statistics, LinearAlgebra, StatsPlots 

include("copem2Daux.jl")

"""
    dstat(v1::Vector, v2::Vector)

Given two vectors of equal length `v1` and `v2` with data that are of (or can be converted to)
`Float` type, a very small random noise is added if there are repeated values, and then a
named tuple with the following keys is generated:

- `summary`: summary marginal and bivariate statistics
- `con`: Empirical Spearman's concordance measure
- `dep`: Empirical Schweizer's dependence measure
- `gini`: Empirical Gini's concordance measure
- `pearson`: Empirical Pearson's correlation
- `xy`: a matrix with the given vectors as columns without repeated values
- `x`, `y`: the given vectors but without repeated values
- `xstat`, `ystat`: named tuples with marginal descriptive statistics of the original vectors
- `Rxy`: 2-column matrix with the marginal ranks of each variable as columns
- `Cxy`: n×n matrix with empirical copula values of n⋅C(i/n, j/n) for i,j ∈ {1,...,n} where n is the sample size
- `fCxy:` 2-argument function that given (i,j) ∈ {1,...,n}^2 returns the value of n⋅C(i/n,j/n) where n is the sample size
- `Dxy`: 3-column matrix where:
> Column 1: values from 0 to n, where n is the sample size

> Column 2: values of n⋅C(i/n,i/n) for i ∈ {0,...,n} (empirical diagonal)

> Column 3: values of n⋅C(i/n, 1-i/n) for i ∈ {0,...,n} (empirical secondary diagonal)


## Example 
```
x = rand(300);
y = @. x*(1-x) + 0.2*rand();
d = dstat(x, y);
keys(d)
d.summary
d.con, d.dep, d.gini, d.pearson
d.Rxy
d.Cxy
d.Cxy[122, 245]
d.fCxy(122, 245)
d.Dxy
d.Cxy[100, 100], d.Cxy[100, 200]
println(d.Dxy[101, :])
```
"""
function dstat(v1::Vector, v2::Vector)
    if length(v1) ≠ length(v2)
        @error "Vectors must be of equal length."
        return nothing
    end
    try
        v1 = float(v1)
        v2 = float(v2)
    catch
        @error "Data cannot be converted to Float type."
        return nothing
    end    
    r1 = resumen(v1)
    r2 = resumen(v2)
    desempate!(v1)
    desempate!(v2)
    matXY = [v1 v2]
    R = rangobs(matXY)
    C = copem(matXY)
    fC(i, j) = fcopem(i, j, R)
    D = diagem(C)
    δ = condepcopem(C)
    γ = condiagem(D)
    p = cor(v1, v2)
    t1 = append!(pushfirst!(collect(keys(r1)), :variable), [:Spearman, :Schweizer, :Gini, :Pearson])
    t2 = vcat(['X'], collect(r1), [δ.con, δ.dep, γ.gini, p])
    t3 = vcat(['Y'], collect(r2), [δ.con, δ.dep, γ.gini, p])
    T = [t1 t2 t3]
    salida = (summary = T, xy = matXY, x = v1, y = v2, xstat = r1, ystat = r2, 
              Rxy = R, Cxy = C, fCxy = fC, Dxy = D, con = δ.con, 
              dep = δ.dep, gini = γ.gini, pearson = cor(v1, v2)
    )
    if r1.repeated > 0 || r2.repeated > 0
        @warn "Variable 1: $(r1.repeated) repeated values jittered out of $(r1.samplesize)"
        @warn "Variable 2: $(r2.repeated) repeated values jittered out of $(r2.samplesize)"
        println()
    end
    return salida 
end


"""

    dplot(d; regfun = false, cregfun = false, w = 500, h = 450)

Given a named tuple `d` generated by the function `dstat` generates a 3×3 matrix 
with the following plots, according to their position in the matrix:

- [2,2] scatterplot of the original bivariate data but without repeated values
- [1,2] scatterplot of the ranks of the bivariate data (without repeated values)
- [1,3] graph of the diagonal of the empirical copula along with copula bounds and independence
- [2,3] graph of the secondary diagonal of the empirical copula along with copula bounds and independence
- [3,1] barplot of the values of the absolute value of Spearman's conocordance and Schweizer's dependence. If the color of Spearman's bar is red means that it's value is negative
- [3,2] absolute frecuency histogram of the first variable
- [3,3] boxplot of the first variable
- [2,1] absolute frecuency histogram of the second variable
- [1,1] boxplot of the second variable variable

The optional named parameter `regfun` (by default is `false`) if `true` then a linear regression curve
is added to the original data scatterplot, or you may provide a regression function `f(a,x)` 
where `0<a<1` (quantile). Similar comment for `cregfun` but for the rankplot.

The optional named paramenters `w` and `h` are respectively the width and height (in pixels) of the matrix
plot which if not specified default to `w = 500` and `h = 450`. 
    
The function returns a named tuple of plot objects. If dp = dplot(d) then:

- dp.all = 3×3 matrix plot with all the plots
- dp.scatter = [2,2]
- dp.rank = [1,2]
- dp.diag = [1,3]
- dp.diag2 = [2,3]
- dp.condep = [3,1]
- dp.hist1 = [3,2]
- dp.box1 = [3,3]
- dp.hist2 = [2,1]
- dp.box2 = [1,1]

## Example
```
x = rand(300);
y = @. x*(1-x) + 0.2*rand();
d = dstat(x, y);
dp = dplot(d)
using StatsPlots (external package)
plot(dp.all)
```
"""
function dplot(d; regfun = false, cregfun = false, w = 500, h = 450)
    # `d` a tuple generated by `dstat`
    # auxiliary formulas for simple linear regression
    n = d.xstat.samplesize
    b(x, y) = (n*x⋅y - sum(x)*sum(y)) / (n*sum(x .^ 2) - sum(x)^2)
    a(x, y) = mean(y) - b(x, y)*mean(x)
    # scatterplot
    g22 = scatter(d.x, d.y, legend = false, ms = 1, mc = :black, 
    xtickfont = 5, ytickfont = 5, xlabel = "X", ylabel = "Y")
    if regfun == true # plot linear regression
        β = b(d.x, d.y)
        α = a(d.x, d.y)
        xx = range(d.xstat.min, d.xstat.max, length = 300)
        yy = α .+ β.*xx
        plot!(xx, yy, color = :red, lw = 1.5)
    elseif regfun == false
        nothing 
    else # must be a regression function `f(a,x)` where 0<a<1 (quantile)
        xx = range(d.xstat.min, d.xstat.max, length = 300)
        yy1 = regfun.(0.1, xx)
        yy2 = regfun.(0.5, xx)
        yy3 = regfun.(0.9, xx)
        plot!(xx, yy1, color = :violet, lw = 1.5)
        plot!(xx, yy2, color = :red, lw = 1.5)
        plot!(xx, yy3, color = :violet, lw = 1.5)
    end
    # rankplot # n = d.xstat.samplesize
    g12 = scatter(d.Rxy[:, 1]/n, d.Rxy[:, 2]/n, legend = false, ms = 1, mc = :black, xlims = (0,1), ylims = (0,1), xtickfont = 5, ytickfont = 5)
    # g12 = scatter(d.Rxy[:, 1]/n, d.Rxy[:, 2]/n, legend = false, ms = 1, mc = :violet, xlims = (0,1), ylims = (0,1), xtickfont = 5, ytickfont = 5)
    u = range(0.0, 1.0, length = 300)
    plot!(u, u, lw = 1.0, color = :gray)
    plot!(u, 1.0 .- u, lw = 1.0, color = :gray)
    if cregfun == true # plot linear regression
        β = b(d.Rxy[:, 1]/n, d.Rxy[:, 2]/n)
        α = a(d.Rxy[:, 1]/n, d.Rxy[:, 2]/n)
        uu = range(0.0, 1.0, length = 300)
        vv = α .+ β.*uu
        plot!(uu, vv, color = :red, lw = 1.5)
    elseif cregfun == false
        nothing 
    else # must be a regression function `cf(a,u)` where 0<a<1 (quantile)
        uu = range(0.0, 1.0, length = 300)
        vv1 = cregfun.(0.1, uu)
        vv2 = cregfun.(0.5, uu)
        vv3 = cregfun.(0.9, uu)
        plot!(uu, vv1, color = :violet, lw = 1.5)
        plot!(uu, vv2, color = :red, lw = 1.5)
        plot!(uu, vv3, color = :violet, lw = 1.5)
    end    
    # histograms and boxplots
    g32 = histogram(d.x, legend = false, xtickfont = 5, ytickfont = 5, color = :steelblue3, alpha = 0.7, lw = 0.2, yaxis = :flip, bins = range(minimum(d.x), maximum(d.x), length = 30))
    # g32 = histogram(d.x, legend = false, xtickfont = 5, ytickfont = 5, color = :cyan3, yaxis = :flip, bins = range(minimum(d.x), maximum(d.x), length = 30))
    # g32 = density!(d.x, trim = true, color = :black, lw = 3, yaxis = :flip, xtickfont = 5, ytickfont = 5)
    g33 = boxplot(d.x, xtickfont = 5, ytickfont = 5, color = :steelblue3, alpha = 0.7, lw = 0.7, permute = (:x, :y), whisker_width = :half, xshowaxis = false, bar_widths = 0.1)
    # g33 = boxplot(d.x, xtickfont = 5, ytickfont = 5, color = :cyan3, permute = (:x, :y))
    g21 = histogram(d.y, legend = false, xtickfont = 5, ytickfont = 5, color = :palegreen3, alpha = 0.7, lw = 0.2, permute = (:x, :y), yaxis = :flip, bins = range(minimum(d.y), maximum(d.y), length = 30))
    # g21 = histogram(d.y, legend = false, xtickfont = 5, ytickfont = 5, color = :orange, permute = (:x, :y), yaxis = :flip, bins = range(minimum(d.y), maximum(d.y), length = 30))
    # g21 = density!(d.y, trim = true, color = :orange, lw = 3, permute = (:x, :y), yaxis = :flip, xtickfont = 5, ytickfont = 5)
    g11 = boxplot(d.y, xtickfont = 5, ytickfont = 5, color = :palegreen3, alpha = 0.7, lw = 0.7, whisker_width = :half, xshowaxis = false)
    # g11 = boxplot(d.y, xtickfont = 5, ytickfont = 5, color = :orange)
    # empirical diagonal
    t = collect(range(0.0, 1.0, length = 1_001));
    dW = @. (2t - 1) * (t > 1/2)
    dP = t .^ 2
    dM = t
    plot(t, dM, legend = false, xtickfont = 5, ytickfont = 5, color = :gray)
    plot!(t, dP, color = :gray)
    plot!(t, dW, color = :gray)
    g13 = scatter!(d.Dxy[:, 1]/n, d.Dxy[:, 2]/n, ms = 0.5)
    # empirical secondary diagonal
    d2W = zeros(length(t))
    d2P = @. t*(1-t)
    d2M = @. t*(t ≤ 0.5) + (1-t)*(t > 0.5)
    plot(t, d2M, legend = false, xtickfont = 5, ytickfont = 5, color = :gray)
    plot!(t, d2P, color = :gray)
    plot!(t, d2W, color = :gray)
    g23 = scatter!(d.Dxy[:, 1]/n, d.Dxy[:, 3]/n, ms = 0.5)
    # Spearman and Schweizer measures
    spcolor = d.con ≥ 0.0 ? :black : :lightgoldenrod1
    g31 = bar(["Spearman", "Schweizer"], [abs(d.con), d.dep], fillcolor = [spcolor, :black],
              fillalpha = [0.7,0.7], xtickfont = 5, ytickfont = 5, ylims = (0,1), lw = 0.7)
    # dplot
    dp = plot(g11, g12, g13, g21, g22, g23, g31, g32, g33, 
              layout = (3,3), legend = false, size = (w, h)
    )
    return (all = dp, scatter = g22, rank = g12, diag = g13, diag2 = g23, condep = g31, hist1 = g32, box1 = g33, hist2 = g21, box2 = g11)
end

end