Skip to content

Conversation

@Santymax98
Copy link
Contributor

Okay, let's forget about the pr I made by mistake.

In this PR, I added the vcov function for the fit API. This ensures that it's the user's decision whether to calculate it or not. I also added a function to print the values ​​for copula measurements (again, the decision is up to the user). Here are some examples.

julia> G = GaussianCopula([1.0 0.3 0.2; 0.3 1.0 0.4; 0.2 0.4 1.0])
GaussianCopula{3, Matrix{Float64}}(Σ = [1.0 0.3 0.2; 0.3 1.0 0.4; 0.2 0.4 1.0]))

julia> U = rand(G, 300)
3×300 Matrix{Float64}:
 0.222071  0.991229  0.343356  0.447073  0.756457  0.909128  0.429314    0.117147  0.23917   0.886017  0.779509  0.239024   0.252678  0.0907745
 0.53644   0.790178  0.593496  0.45622   0.763848  0.166268  0.464444     0.51989   0.601099  0.965531  0.710267  0.0481338  0.156809  0.477384
 0.456802  0.832535  0.644807  0.571332  0.370216  0.255189  0.642993     0.57973   0.637748  0.894761  0.607114  0.164993   0.290189  0.67479

julia> fit(CopulaModel, GaussianCopula, U; vcov = true)
[ Info: Running the itau/irho/ibeta routine from the generic implementation
Gaussian d=3 fitted via itau
Number of observations:       300
Null Loglikelihood:        0.0000
Loglikelihood:            32.2020
LR Test (vs indep. copula): 64.40 ~ χ²(3)  =>  p = 6.727e-14
AIC: -58.404       BIC: -47.293
Converged: true   Iterations: 52   Elapsed: 3.077s
────────────────────────────────────────────────────────────────────────────────────────
Parameter          Estimate      Std.Err   z-value   Pr(>|z|)       95% Lo       95% Hi
────────────────────────────────────────────────────────────────────────────────────────
Σ_1_2              0.219403   0.00354438    61.902          0     0.212456      0.22635
Σ_1_3              0.160347   0.00349548    45.873          0     0.153496     0.167198
Σ_2_3               0.38445   0.00319102   120.479          0     0.378196     0.390704
────────────────────────────────────────────────────────────────────────────────────────

For Archimedean Copulas

julia> C = ClaytonCopula(3, 3.5)
ClaytonCopula{3, Float64}(θ = 3.5,)

julia> U = rand(C, 300)
3×300 Matrix{Float64}:
 0.994932  0.476203  0.52375   0.00110445   0.270389  0.579527  0.959464    0.896951  0.858705  0.223188  0.261711  0.924947  0.322537  0.552703
 0.823358  0.389963  0.492272  0.000697532  0.313849  0.905658  0.565444     0.837394  0.81539   0.198364  0.309672  0.46887   0.358373  0.705274
 0.842299  0.945494  0.534202  0.00101964   0.234375  0.495992  0.665326     0.53361   0.866946  0.266491  0.939297  0.57424   0.369273  0.742542

julia> fit(CopulaModel, ClaytonCopula, U; method=:itau, vcov = true)
Archimedean d=3 fitted via itau
Number of observations:       300
Null Loglikelihood:        0.0000
Loglikelihood:           531.8877
LR Test (vs indep. copula): 1063.78 ~ χ²(1)  =>  p = 2.467e-233
AIC: -1061.775       BIC: -1058.072
Converged: true   Iterations: 0   Elapsed: 2.401s
────────────────────────────────────────────────────────────────────────────────────────
Parameter          Estimate      Std.Err      z-value     Pr(>|z|)       95% Lo       95% Hi
────────────────────────────────────────────────────────────────────────────────────────
θ                    3.7830       1.3000       2.9101       0.0036       1.2352       6.3309
────────────────────────────────────────────────────────────────────────────────────────

Copulas with more than 2 params... with fallback to avoid errors

julia> B = BB1Copula(2, 1.5, 3.2)
BB1Copula{2, Float64}(θ = 1.5, δ = 3.2)

julia> U = rand(B, 300)
2×300 Matrix{Float64}:
 0.0868505  0.369871  0.0876678  0.192073  0.902674  0.411426  0.0638137    0.742172  0.583333  0.532168  0.0187803  0.0129811  0.554113  0.950487
 0.150734   0.595093  0.139971   0.145485  0.892565  0.263838  0.105728      0.884073  0.539801  0.551795  0.010793   0.0161772  0.713115  0.93326

julia> fit(CopulaModel, BB1Copula, U; vcov = true)
┌ Warning: vcov(hessian) failed (NaN/Inf). Falling back to jackknife.
└ @ Copulas ~/Desktop/Repositorios/Copulas.jl/src/Fitting.jl:490
Archimedean d=2 fitted via mle
Number of observations:       300
Null Loglikelihood:        0.0000
Loglikelihood:           409.2012
LR Test (vs indep. copula): 818.40 ~ χ²(2)  =>  p = 1.933e-178
AIC: -814.402       BIC: -806.995
Converged: true   Iterations: 9   Elapsed: 3.916s
────────────────────────────────────────────────────────────────────────────────────────
Parameter          Estimate      Std.Err      z-value     Pr(>|z|)       95% Lo       95% Hi
────────────────────────────────────────────────────────────────────────────────────────
θ                    1.2132       0.0117     103.3044       0.0000       1.1902       1.2362
δ                    3.5206       0.0154     228.2639       0.0000       3.4903       3.5508
────────────────────────────────────────────────────────────────────────────────────────

For SklarDist

julia> X₁, X₂, X₃ = Gamma(2,3), Pareto(), LogNormal(0,1)
(Gamma{Float64}(α=2.0, θ=3.0), Pareto{Float64}(α=1.0, θ=1.0), LogNormal{Float64}(μ=0.0, σ=1.0))

julia> C = ClaytonCopula(3, 1.7)
ClaytonCopula{3, Float64}(θ = 1.7,)

julia> D = SklarDist(C,(X₁,X₂,X₃))
SklarDist{ClaytonCopula{3, Float64}, Tuple{Gamma{Float64}, Pareto{Float64}, LogNormal{Float64}}}(
C: ClaytonCopula{3, Float64}(θ = 1.7,)
m: (Gamma{Float64}(α=2.0, θ=3.0), Pareto{Float64}(α=1.0, θ=1.0), LogNormal{Float64}(μ=0.0, σ=1.0))
)


julia> x = rand(D,1000)
3×1000 Matrix{Float64}:
 2.05448  2.3674   13.5994    4.19579  13.6778   1.47767   2.26288   1.96553     7.10346   7.15242  12.8446  2.50478  6.31336   6.18693  7.3298
 1.21988  1.55024   4.87178   1.84246   1.27558  1.02665   1.41345   1.69677      1.42585   4.57979   2.8184  1.88363  3.5216    2.29623  1.2645
 3.6093   2.01363   0.702068  1.43357   1.13343  0.788215  0.409093  0.613827     0.544334  1.70939   1.2616  1.52514  0.653832  1.03131  0.334373

julia> fit(CopulaModel,
           SklarDist{ClaytonCopula,Tuple{Gamma,Pareto,LogNormal}},
           x; copula_method=:mle, sklar_method=:ifm,
           copula_kwargs=(vcov=true,), derived_measures=true)
SklarDist{Copula=Archimedean d=3, Margins=(Gamma, Pareto, LogNormal)} fitted via mle
Number of observations:      1000
Null Loglikelihood:    -6176.8206
Loglikelihood:         -5423.4946
LR Test (vs indep. copula): 1506.65 ~ χ²(1)  =>  p = 0
AIC: 10860.989       BIC: 10895.343
Converged: true   Iterations: 4   Elapsed: 3.540s
──────────────────────────────────────────────────────────
[ Copula ]
──────────────────────────────────────────────────────────
Family: Archimedean d=3
────────────────────────────────────────────────────────────────────────────────────────
Param            Estimate      Std.Err   z-value   Pr(>|z|)       95% Lo       95% Hi
────────────────────────────────────────────────────────────────────────────────────────
θ                  1.5829       0.0075   212.179          0       1.5683       1.5975
────────────────────────────────────────────────────────────────────────────────────────
[ Copula Derived measures ]
Kendall τ(θ)   = 0.4418
Spearman ρ(θ)  = 0.6149
Blomqvist β(θ) = 0.1550
Gini γ(θ)      = 0.5124
Upper λᵤ(θ)    = 0.0000
Lower λₗ(θ)    = 0.4995
Entropy ι(θ)   = -0.7601
──────────────────────────────────────────────────────────
[ Marginals ]
──────────────────────────────────────────────────────────
Margin Dist         Param       Estimate      Std.Err       95% CI
#1     Gamma        α             1.9651            —            —
       Gamma        θ             3.0329            —            —
#2     Pareto       α             0.9651            —            —
       Pareto       θ             1.0003            —            —
#3     LogNormal    μ             0.0115            —            —
       LogNormal    σ             1.0008            —            —

In this last case, I had trouble getting the marginal vcov...

@lrnv
Copy link
Owner

lrnv commented Oct 2, 2025

This looks great, although I think it would be better (more generic) if everything happened in the outer level, without having to go inside the different methods, so that we get vcov whatever the fitting method. I will take a look tomorrow to see what I can do.

@codecov
Copy link

codecov bot commented Oct 2, 2025

Codecov Report

❌ Patch coverage is 70.33639% with 194 lines in your changes missing coverage. Please review.
✅ Project coverage is 78.55%. Comparing base (29108e3) to head (fc12c12).
⚠️ Report is 3 commits behind head on main.

Files with missing lines Patch % Lines
src/show.jl 0.00% 129 Missing ⚠️
src/Fitting.jl 84.00% 32 Missing ⚠️
src/utils.jl 58.82% 7 Missing ⚠️
src/WilliamsonTransforms.jl 84.61% 6 Missing ⚠️
src/Generator.jl 86.20% 4 Missing ⚠️
src/MiscellaneousCopulas/FGMCopula.jl 91.11% 4 Missing ⚠️
src/ExtremeValueCopula.jl 83.33% 2 Missing ⚠️
src/Generator/BB2Generator.jl 88.88% 2 Missing ⚠️
src/MiscellaneousCopulas/IndependentCopula.jl 0.00% 2 Missing ⚠️
src/MiscellaneousCopulas/MCopula.jl 0.00% 2 Missing ⚠️
... and 4 more
Additional details and impacted files
@@            Coverage Diff             @@
##             main     #311      +/-   ##
==========================================
+ Coverage   74.50%   78.55%   +4.04%     
==========================================
  Files          84       85       +1     
  Lines        4837     4994     +157     
==========================================
+ Hits         3604     3923     +319     
+ Misses       1233     1071     -162     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@Santymax98
Copy link
Contributor Author

I understand, but sometimes one technique is used for MLE and another for moment methods, finally the safe fallback to light Jackknife because I think a bootstrap would be too heavy

@lrnv
Copy link
Owner

lrnv commented Oct 3, 2025

So if one technique is for MLE and the other one for moment methods, then we could simply dipatch on the method, or even add a new argument vcov_method to decide what we'll do, whihc could have a smart default behavior:

  • If the method is MLE, dispatch by defautl to the MLE hessian
  • If the method is moment-matching, dispatch to the correct code
  • Else dispatch to jacknife by default.

Let me look at your code and try to do that.

Copy link
Owner

@lrnv lrnv left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Santymax So I think i have it working now :

julia> U = rand(GaussianCopula([1.0 0.3 0.2; 0.3 1.0 0.4; 0.2 0.4 1.0]), 300);
julia> fit(CopulaModel, GaussianCopula, U)
────────────────────────────────────────────────────────────────────────────────
[ CopulaModel: Gaussian d=3 ]
────────────────────────────────────────────────────────────────────────────────
Method:                mle
Number of observations: 300
────────────────────────────────────────────────────────────────────────────────
[ Fit metrics ]
────────────────────────────────────────────────────────────────────────────────
Null Loglikelihood:          0.0000
Loglikelihood:              50.4789
LR (vs indep.):        100.96 ~ χ²(3)    p = <1e-16
AIC:                   -94.958
BIC:                   -83.847
Converged:             true
Iterations:            0
Elapsed:               0.000s
────────────────────────────────────────────────────────────────────────────────
[ Parameters ]
────────────────────────────────────────────────────────────────────────────────
Parameter    Estimate   Std.Err   z-value    p-val     95% Lo     95% Hi
Σ₁₂            0.2793    0.0538     5.196 2.038e-07     0.1740     0.3847
Σ₁₃            0.2057    0.0569     3.615 0.0003009     0.0942     0.3173
Σ₂₃            0.4754    0.0435    10.924   <1e-16     0.3901     0.5607


julia>
ulia> U = rand(ClaytonCopula(3, 3.5), 300);
julia> fit(CopulaModel, ClaytonCopula, U; method=:itau)
────────────────────────────────────────────────────────────────────────────────
[ CopulaModel: Archimedean d=3 ]
────────────────────────────────────────────────────────────────────────────────
Method:                itau
Number of observations: 300
────────────────────────────────────────────────────────────────────────────────
[ Fit metrics ]
────────────────────────────────────────────────────────────────────────────────
Null Loglikelihood:          0.0000
Loglikelihood:             456.7272
LR (vs indep.):        913.45 ~ χ²(1)    p = <1e-16
AIC:                   -911.454
BIC:                   -907.751
Converged:             true
Iterations:            0
Elapsed:               1.165s
────────────────────────────────────────────────────────────────────────────────
[ Parameters ]
────────────────────────────────────────────────────────────────────────────────
Parameter    Estimate   Std.Err   z-value    p-val     95% Lo     95% Hi
θ              3.5488    0.0000 3548815.824   <1e-16     3.5488     3.5488


julia>
julia> U = rand(BB1Copula(2, 1.5, 3.2), 300);
julia> fit(CopulaModel, BB1Copula, U)
────────────────────────────────────────────────────────────────────────────────
[ CopulaModel: Archimedean d=2 ]
────────────────────────────────────────────────────────────────────────────────
Method:                mle
Number of observations: 300
────────────────────────────────────────────────────────────────────────────────
[ Fit metrics ]
────────────────────────────────────────────────────────────────────────────────
Null Loglikelihood:          0.0000
Loglikelihood:             376.9690
LR (vs indep.):        753.94 ~ χ²(2)    p = <1e-16
AIC:                   -749.938
BIC:                   -742.530
Converged:             true
Iterations:            8
Elapsed:               2.469s
────────────────────────────────────────────────────────────────────────────────
[ Parameters ]
────────────────────────────────────────────────────────────────────────────────
Parameter    Estimate   Std.Err   z-value    p-val     95% Lo     95% Hi
θ              1.5194    0.2542     5.978 2.259e-09     1.0213     2.0176
δ              2.9502    0.2617    11.275   <1e-16     2.4373     3.4630


julia>
julia> x = rand(SklarDist(ClaytonCopula(3, 1.7),(Gamma(2,3), Pareto(), LogNormal(0,1))),1000);
julia> fit(CopulaModel, SklarDist{ClaytonCopula,Tuple{Gamma,Pareto,LogNormal}}, x)
────────────────────────────────────────────────────────────────────────────────
[ CopulaModel: SklarDist (Copula=Archimedean d=3, Margins=(Gamma, Pareto, LogNormal)) ]
────────────────────────────────────────────────────────────────────────────────
Copula:                Archimedean d=3
Margins:               (Gamma, Pareto, LogNormal)
Methods:               copula=mle, sklar=ifm
Number of observations: 1000
────────────────────────────────────────────────────────────────────────────────
[ Fit metrics ]
────────────────────────────────────────────────────────────────────────────────
Null Loglikelihood:      -6062.1153
Loglikelihood:                 -Inf
LR (vs indep.):        -Inf ~ χ²(1)    p = 1
AIC:                   Inf
BIC:                   Inf
Converged:             true
Iterations:            1
Elapsed:               3.258s
────────────────────────────────────────────────────────────────────────────────
[ Dependence metrics ]
────────────────────────────────────────────────────────────────────────────────
Kendall τ:             0.4290
Spearman ρ:            0.5995
Blomqvist β:           0.1464
────────────────────────────────────────────────────────────────────────────────
[ Copula parameters ] (vcov=hessian)
────────────────────────────────────────────────────────────────────────────────
Parameter    Estimate   Std.Err   z-value    p-val     95% Lo     95% Hi
θ              1.5026    0.0515    29.184   <1e-16     1.4017     1.6035
────────────────────────────────────────────────────────────────────────────────
[ Marginals ]
────────────────────────────────────────────────────────────────────────────────
Margin Dist       Param    Estimate   Std.Err 95% CI
#1     Gamma      α          2.1972    0.0918 [2.0173, 2.3770]
                  θ          2.7371    0.1284 [2.4855, 2.9887]
#2     Pareto     α          0.9967    0.5524 [-0.0859, 2.0794]
                  θ          1.0016    0.5560 [-0.0881, 2.0914]
#3     LogNormal  μ         -0.0028    0.0313 [-0.0641, 0.0585]
                  σ          0.9886    0.0221 [0.9452, 1.0319]


julia> 

Now by default everything is computed, and i have tried to make the output better-looking. Also I have fixed the marginal vcov.

However this require a few tests: there are now many possibilities (fitting method for the copula, for the sklar, vcov method, etc...) and thus we need to test all of them. I will now review once more the code to try to simplify it again, if you could propose :

1° tests that corresponds to these new functionalities
2° Proper documentation of what the show() methods outputs

That would be nice :)

@lrnv
Copy link
Owner

lrnv commented Oct 7, 2025

@Santymax98 Okay this is correctly merged now ! Do not forget to ]resolve because the packages dependencies changed. I will look at the tests now.

@lrnv
Copy link
Owner

lrnv commented Oct 7, 2025

@Santymax98 : Regarding specific methods for the BBCops generator derivatives, I have simply commented out your code for the moment, since it is not the main question of this PR. Once this one is merged, if you want to test this idea again please in another pull-request. This is already a very messy pull-request no need to add more.

Let's see what the tests say

@Santymax98
Copy link
Contributor Author

I completely agree. It's not essential, and perhaps we can run some benchmarks to see if it's better or not at the appropriate time. For my part, everything's fine.

@lrnv
Copy link
Owner

lrnv commented Oct 8, 2025

Okay. In the mean time, i would like to validate the dependence metrics to be certain that they work on every model we have. I will therefore push a nex tests into the generics and we'll see :)

It will probably fail on a lot of cases, but thats good we have to fix them

@Santymax98
Copy link
Contributor Author

I think the same and maybe polish the tail estimators and theorical tails...

@lrnv
Copy link
Owner

lrnv commented Oct 8, 2025

@Santymax98 Okay I think we are good to merge now. The only missing thing is correct documentation of the fitting interface and the vcov things. I was going to do it, but if you could do it that would be nice.

I will rather concentrate myself on 1) producing benchmarks so that we can see where are the problematic points and 2) make the testing framework a bit faster, 45minutes is painfully slow...

@lrnv lrnv merged commit 0463930 into lrnv:main Oct 8, 2025
5 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants