Stata and Mata have no expm1(x) function, i.e. a function which computes the function exp(x)-1.

Computing exp(x)-1 suffers from cancellation errors for small values of x and will be inaccurate in this case (see in particular this paper https://cran.r-project.org/web/packa...1mexp-note.pdf).

This function is in particular useful to compute accurately 1-exp(-x) which looks like the invcloglog functon (1-exp(-exp(x)). We can compute 1-exp(-x)=-expm1(-x).

It also seems that the invcloglog implementation in Stata and Mata does not use this function. The mata code for invcloglog(x) basically returns 1-exp(-exp(x)) for all values, only for non-missing values of x>500 the result is set equal to 1.

Note that the function expm1() is part of the C and C++ language as well as the function log1p().

I have tried to implement the function myself in Mata but I am in doubt whether I am doing it right. Basically my approach is to use a Taylor expansion of order two for values less or equal to 1e-5. I am not sure whether computing more terms will actually help.

I have read somewhere else that one can use exp(x/2)*sinh(x/2) for values of x between epsilon()/2 and log(2) (where epsilon is the machine precision) and for values less than epsilon/2 then expm1 = x.

But I don't know which approach is right, if there is any.

I am no expert in this matter and I would very much like to benefit from your insights about this issue.

Does someone know which approximation I should use (as well as the thresholds ) or could give som advice or inputs in how to compute this function?

Best

Christophe

Computing exp(x)-1 suffers from cancellation errors for small values of x and will be inaccurate in this case (see in particular this paper https://cran.r-project.org/web/packa...1mexp-note.pdf).

This function is in particular useful to compute accurately 1-exp(-x) which looks like the invcloglog functon (1-exp(-exp(x)). We can compute 1-exp(-x)=-expm1(-x).

It also seems that the invcloglog implementation in Stata and Mata does not use this function. The mata code for invcloglog(x) basically returns 1-exp(-exp(x)) for all values, only for non-missing values of x>500 the result is set equal to 1.

Note that the function expm1() is part of the C and C++ language as well as the function log1p().

I have tried to implement the function myself in Mata but I am in doubt whether I am doing it right. Basically my approach is to use a Taylor expansion of order two for values less or equal to 1e-5. I am not sure whether computing more terms will actually help.

Code:

real matrix expm1(real matrix x) { real matrix index index = abs(x) :<= 1e-5 return(!index:*(exp(x):-1):+index:*(x+0.5*x:^2)) }

I have read somewhere else that one can use exp(x/2)*sinh(x/2) for values of x between epsilon()/2 and log(2) (where epsilon is the machine precision) and for values less than epsilon/2 then expm1 = x.

But I don't know which approach is right, if there is any.

I am no expert in this matter and I would very much like to benefit from your insights about this issue.

Does someone know which approximation I should use (as well as the thresholds ) or could give som advice or inputs in how to compute this function?

Best

Christophe

## Comment