  At present, there is still no officially accepted and extensively verified
implementation of computing the gamma difference distribution allowing unequal
shape parameters. We explore four computational ways of the gamma difference
distribution with the different shape parameters resulting from time series
kriging, a forecasting approach based on the best linear unbiased prediction,
and linear mixed models. The results of our numerical study, with emphasis on
using open data science tools, demonstrate that our open tool implemented in
high-performance Python(with Numba) is exponentially fast, highly accurate, and
very reliable. It combines numerical inversion of the characteristic function
and the trapezoidal rule with the double exponential oscillatory transformation
(DE quadrature). At the double 53-bit precision, our tool outperformed the
speed of the analytical computation based on Tricomi's $U(a, b, z)$ function in
CAS software (commercial Mathematica, open SageMath) by 1.5-2 orders. At the
precision of scientific numerical computational tools, it exceeded open SciPy,
NumPy, and commercial MATLAB 5-10 times. The potential future application of
our tool for a mixture of characteristic functions could open new possibilities
for fast data analysis based on exact probability distributions in areas like
multidimensional statistics, measurement uncertainty analysis in metrology as
well as in financial mathematics and risk analysis.
