BiCopTau2Par for Frank

[oezgesahin]

Hello, In the estimation of the bivariate Frank copula parameter via Kendall's tau inversion, tau values being larger than 0.96 is causing an error, e.g.,:

BiCopTau2Par(5, 0.966)
Error in uniroot(function(x) tau - frankTau(x), lower = 0 + .Machine$double.eps^0.5,:f() values at end points not of opposite sign

The problem seems to result from the values frankParGrid and frankTauVals in BiCopPar2Tau that should be edited according to #74?

[tnagler]

I think you're right. Could you submit a pull request with a fix?

[oezgesahin]

Sure, but I am unsure from a theoretical point of view. For the bivariate Frank copula with parameter theta, it holds: tau = 1 - 4/theta + 4/theta*D(theta), where D: debye function. So, for theta=500, I am obtaining tau=0.99:

frankPar <- 500
tau <- 1 - 4/frankPar + 4/frankPar * copula::debye1(frankPar)
tau 
[1] 0.9920263

It implies that the maximum tau value with our upper bound (35) would be around 0.89.

frankPar <- 35
tau <- 1 - 4/frankPar + 4/frankPar * copula::debye1(frankPar)
tau 
[1] 0.8910855

How should we proceed?

[tnagler]

Actually, the easiest should be to restrict the parameter range when root solving here. Right now, a parameter > 35 could result from calcPar(), but it would be truncated immediately after here. So truncating already in calcPar() wouldn't make things worse. In that case, we could also recompute the lookup table with a narrower range.

What do you think?

[oezgesahin]

Thanks! I think that we could make the grid and its corresponding values wider so that it can compute the parameter-tau relationship until tau=0.99 with the theoretical formula. Accordingly, we should extend the parameter range in root solving to make sure that lower and upper bounds return opposite signs as required. For any tau value > 0.89, the returned parameter value will be 35 in this scenario. Does that make sense?

[oezgesahin]

Also, do you know why these (1, 2) are commented and values are written one by one instead?

[tnagler]

But what's the point in making the grid wider, if the result is 35 anyway? Can't we just restrict the root solving range then?

We store the values one by one, so VineCopula does not depend on the copula package (and as a result also gsl, which requires extra system libraries on linux and mac).

[oezgesahin]

If I am not mistaken: approx is interpolating given the grid values for par (x) and tau (y). Therefore, max tau value it can work is bounded by the upper limit of the tau grid. Thus, when we restrict root solving range without making the grid wider, uniroot function could still return positive (tau-interpolation) values for tau > 0.96 in the upper bound, which prevents it from finding a root. Making the grids wider up to tau = 0.99 and adjust the upper bound (app. 500) accordingly should approximate the theoretical relationship with one difference that any parameter > 35 will be truncated to 35. Without adjusting the grid but changing uniroot range:

frankParGrid <- seq(-100, 100, l = 100)
frankTauVals <- 1 - 4/frankParGrid + 4/frankParGrid * copula::debye1(frankParGrid)
tau <- 0.99
v <- uniroot(function(x) tau - frankTau(x),
              lower = 0 + .Machine$double.eps^0.5, upper = 35,
              tol = .Machine$double.eps^0.5)$root
Error in uniroot(function(x) tau - frankTau(x), lower = 0 + .Machine$double.eps^0.5,  : 
  f() values at end points not of opposite sign

Suggestion:

frankParGrid <- seq(-500, 500, l = 100)
frankTauVals <- 1 - 4/frankParGrid + 4/frankParGrid * copula::debye1(frankParGrid)
tau <- 0.99
v <- uniroot(function(x) tau - frankTau(x),
             lower = 0 + .Machine$double.eps^0.5, upper = 500,
             tol = .Machine$double.eps^0.5)$root
# v will be truncated to 35
v
[1] 398.3634

[tnagler]

Ah yeah, I get it. Then even better:

frankParGrid <- c(-10^(5:2), seq(-36, 36, l = 100), 10^(2:5))
frankTauVals <- 1 - 4/frankParGrid + 4/frankParGrid * copula::debye1(frankParGrid)
frankTau <- function (par) approx(x = frankParGrid, y = frankTauVals, xout = par)$y
tau <- 0.99
v <- uniroot(function(x) tau - frankTau(x),
             lower = 0 + .Machine$double.eps^0.5, upper = 1e5,
             tol = .Machine$double.eps^0.5)$root
# v will be truncated to 35
v

That way we don't lose accuracy in the [-35,35] range but can handle taus up to 0.9999.

[oezgesahin]

Agreed, great!

[tnagler]

solved by #79