PDF returning Inf when should be zero with Frank Copula

[rand5]

MWE:

θ = -2.5
d= 2

C = FrankCopula(d,θ)

D1 = SklarDist(C,(Normal(-2.,1),Normal(-0.3,0.1)))

pdf(D1, [2.,-2.]) # = Inf

Returns 0 if θ = 2.5 instead of θ = -2.5

[lrnv]

Are you sure it should be zero ? I guess that inverting the parameter basically rotates the frank copula by 1/4 turn, so it appears logic that the limit in the lower corner is not the same. Do you have a formal reference for this pdf ?

[lrnv]

From this wikipedia page, after checking that the Frank copula generator is the same as ours, we get for the denisty the formula

$$c(u,v) = \frac{-\theta e^{-\theta(u+v)}(e^{-\theta}-1)}{(e^{-\theta}-e^{-\theta u}-e^{-\theta v}+e^{-\theta(u+v)})^2}$$

and thus

$$c(0,0) = \frac{-\theta (e^{-\theta}-1)}{(e^{-\theta}-1)^2} = \frac{\theta}{1-e^{-\theta}}$$

which is increasing in theta, one if theta=0, and tends to 0 at theta = -infty and +infty as theta -> +infty. So basically it says that the value of the pdf in zero should not be infinity but rather this ? Maybe there is a mistake somewhere

[lrnv]

Moreover, the R::copula package seems to return 0 in both cases, while the limiting behavior does not look like 0:

> dCopula(c(0,0), frankCopula(2.5))
[1] 0
> dCopula(c(0,0), frankCopula(-2.5))
[1] 0
> dCopula(c(0.00001,0), frankCopula(-2.5))
[1] 0
> dCopula(c(0.00001,0.00001), frankCopula(-2.5))
[1] 0.2235749
> dCopula(c(0.000001,0.000001), frankCopula(-2.5))
[1] 0.2235648
> dCopula(c(0.000001,0.000001), frankCopula(2.5))
[1] 2.72355
> 

So I think we have two options: Either we keep the current behavior, or we fix it by using the $\frac{\theta}{1 - e^{-\theta}}$, which seem to correspond to the values given by R near (0,0):

julia> θ = 2.5
2.5

julia> θ / (1-exp(-θ))
2.72356372458463

julia> θ = -2.5
-2.5

julia> θ / (1-exp(-θ))
0.22356372458463003

julia>

[lrnv]

We could also fix the other corners for $c(0,1), c(1,0)$ and $c(1,1)$ from the formula. What would be interesting is to try to fix it automatically also for higher-dimensional frank copulas... and that i really dont knw how to do it. But we can still use this formula as a test case (or even use it directly by writting a bivariate density method as we did for the Gumbel cdf).

[rand5]

I'm having trouble seeing that this is an issue with the corners. What struck me was not that the results were different if θ = -2.5 versus θ = 2.5, but that the result for θ = -2.5 at that particular location was Inf. I probably don't have the right intuition, but I was thinking if the copula pdf is a valid, the total area under it should be 1, thus at no point should it be Inf.

[lrnv]

I probably don't have the right intuition, but I was thinking if the copula pdf is a valid, the total area under it should be 1, thus at no point should it be Inf.

Simply think about the pdf of a exponential random variable, which is inf at 0: being inf is not a problem for the integral of the function to be 1.

[rand5]

That is a good point.

I just don't see any reason for it to blow up to Inf in the region (approx 1.14 to 4.15) it does based on the equation:

θ = -2.5
d= 2

C = FrankCopula(d,θ)

D1 = SklarDist(C,(Normal(-2.,1),Normal(-0.3,0.1)))

xs = collect(range(-10.,10.,1000))
y = map(xs) do x
    pdf(D1, [x,-x])
end
plot(scatter(;x=xs,y=y))

[lrnv]

I do not even understand what is this graph supposed to be ?

[rand5]

Sorry for the confusion--I'm taking the equation for the convolution, $f_{Z}(z) = \int_{-\infty}^{\infty} f_{XY}(x,z-x)dx$, and looking at the integrand when setting $z=0$, i.e. $f_{XY}(x,-x)$ as a function of $x$. If one were to try to compute the convolution at $z=0$, the integral doesn't work because there is this region of $f_{XY}(x,-x)$ where the pdf (pdf(D1, [x,-x])) is returning Inf between approx x= 1.14 to 4.15.

[lrnv]

Hey @rand5 I think that all this should be fixed in the last version (0.1.31). Can you check again to tell me if you still have issues ?