Fourier inverting a characteristic function when we know the resulting pdf has a dirac delta at zero

Question

I have a characteristic function which is associated with a kind of compound Poisson process. $$\phi(u)=\left(\frac{k - iu}{k - iu\mathrm{e}^{-\eta(t-t_0)}}\right)^{-\lambda/\eta}$$ I want to numerically calculate its probability density function, but I am having trouble due to the singularity at zero, which vanishes as $t \rightarrow +\infty$, as the distribution tends to the gamma distribution $\Gamma(\lambda/\kappa, \eta)$. So when $t>>t_0$ I am already using this approximation and it has good results. Now I need to fix this for the non-limit case.

So I want to obtain

$$f(x) = \frac{1}{2\pi}\int_{-\infty}^{+\infty}\mathrm{e}^{-iux}\phi(u)du$$

My current idea is to try to find a way to obtain a characteristic function $\phi_g(u)$ that is a version of $\phi(u)$ with the Dirac delta contribution removed. Afterwards, I would numerically invert $\phi_g(u)$ and perform the appropriate adjustments when integrating with respect to $f(x)$ for calculating probability convolutions.

My naive attempt

I know that the singularity at zero comes from the probability of the compound Poisson process having no jumps, which is the probability of a Poisson random variable with parameter $\lambda(t-t_0)$ being zero:

$$P[N=0]=\mathrm{e}^{-\lambda(t-t_0)}$$

I think another way of seeing this might be through the limit

$$\mu_f(\{0\})=\displaystyle \lim_{R\rightarrow +\infty} \frac{1}{2\pi R} \int_{-R}^{R}\mathrm{e}^{-iu\cdot 0}\phi(u)du$$

Which I am currently trying to solve, and would expect it to give me $\mathrm{e}^{-\lambda(t-t_0)}$.

I did a substitution $k-iu=y$ and got to

$$\mathrm{e}^{-\lambda(t-t_0)}\int_{Y_{-}}^{Y_{+}}i\left(\frac{-y-ik}{y}\right)^{\lambda/\eta}dy$$

And I've been trying a lot of things but I can't seem to solve this. Although I'm not sure if there's a way of getting the limit of this integral without having to find its expression.

I believe this would mean that $f(x)$ would be

$$f(x) = \mathrm{e}^{-\lambda(t-t_0)}\delta (x) + g(x)$$

And now I'm going to abuse this by using the characteristic function of the Dirac delta and get

$$\phi(u) = \int_{-\infty}^{+\infty}\mathrm{e}^{iux}f(x)dx$$

$$\phi(u) = \int_{-\infty}^{+\infty}\mathrm{e}^{iux}\left(\mathrm{e}^{-\lambda(t-t_0)}\delta(x) + g(x)\right)dx$$

$$\phi(u) = \mathrm{e}^{-\lambda(t-t_0)} + \int_{-\infty}^{+\infty}\mathrm{e}^{iux}g(x)dx$$

$$\phi_g(u) = \phi(u) - \mathrm{e}^{-\lambda(t-t_0)}= \left(\frac{k - iu}{k - iu\mathrm{e}^{-\eta(t-t_0)}}\right)^{-\lambda/\eta} - \mathrm{e}^{-\lambda(t-t_0)}$$ And at the end I want to add both contributions (the $\delta$ and the inverted $\phi_g$) to obtain the full pdf $f(x)$. Is this approach valid? When I try to invert this $\phi_g(u)$ numerically, I still get a spike at zero. I'll try to post some plots later.

My adequacy test for this would be to get the integral of the pdf $f(x)$ I obtained to be 1, but it's far from it.

Using complex analysis?

Now, looking at the original characteristic function, and ignoring this information about expecting the singularity to be at zero, I get confused to due my lack of sufficient knowledge on complex analysis.

It seems we have a pole at $u = -ik$, which is removable in the $t \rightarrow +\infty$ limit. However, I am not sure how this could be used for my purposes. I'm not confident we can actually analytically calculate this full pdf by calculating the residues around this pole, but I'm not sure if there's a way to just remove the contribution of the Dirac delta at zero and get an expression for the rest (which I could try to invert numerically).

Thanks and sorry if the question is badly posed due to my poor knowledge on complex analysis.

Related question that might help give some background(Characteristic function of Gamma-OU process (cross-posted from quant.stackexchange))

Numerical Tests

I have made a few tests that seem to indicated the approach of subtracting the $\mathrm{e}^{-\lambda(t-t_0)}\delta(x)$ component to the characteristic function $\phi(u)$ seems to work well for the cases where the spike is significant. As mentioned below in a comment by @sb945, I thought the spikes persisting could have been caused by $\lambda/\eta =1$, but it seems to work for that case as well (surprisingly). What seemed to be causing some numerical instability was the fact that initially (to avoid the zero) I had numerically inverted the characteristic function for an interval starting an machine epsilon in front of zero, to avoid it. The correction doesn't seem to be applied unless we invert $\phi_g(u)$ starting exactly on zero. I have inverted the characteristic function using a quadrature/fft algorithm using 2^15 points, I choose the interval to be between zero and the mean plus 10 standard deviations of the random variable. I obtain an array with the values of $\phi_g(u)$ that I interpolate linearly to obtain the numerical proxy for it. I check if it's normalized by doing numerical integration (gauss kronod quadrature) and adding the mass at zero. If anyone is interested, I can share the code (it's in Julia).

First test case

$k=1; \lambda=1; \eta=1; t_0=0; t=1$

probability of zero = 0.36787944117144233

integral of inverted $\phi_g$ = 0.6320861572332717

integral of inverted $\phi$ = 0.8160258787203427

full probability (using $\phi_g$ and $\delta$) = 0.999965598404714

In this case with $\lambda/\eta=1$ it performs well.

Second test case

$k=1; \lambda=0.5; \eta=1; t_0=0; t=10$

probability of zero = 0.006737946999085467

integral of inverted $\phi_g$ = 0.9967154997905739

integral of inverted $\phi$ = 1.0000844732901182

full probability (using $\phi_g$ and $\delta$) = 1.0034534467896594

We see that in this case where the probability of zero is small and $\lambda/\eta<1$, inverting the original $\phi$ seems to respect the normalization more.

Third test case

$k=1; \lambda=1; \eta=0.5; t_0=0; t=10$

probability of zero = 4.5399929762484854e-5

integral of inverted $\phi_g$ = 0.9999519102799828

integral of inverted $\phi$ = 0.9999746102448701

full probability (using $\phi_g$ and $\delta$) = 0.9999973102097452

In this case where the probability of zero is also small but $\lambda/\eta>1$, my approach also normalizes close to 1.

Conclusion (for now...)

It seems this approach can give us some decent results numerically. However, I am not going to close this question by answering it myself because I haven't been able to prove that the probability of zero jumps corresponds to that limit I wrote above.

Answer 1

The value of the Fourier transform at 0 is the value of the integral of the original function. It looks like the integrability of your $\phi_g$ will still depend on $\lambda/\eta$ even if $e^{-\lambda(t-t_0)}$ were the correct constant to be subtracted from $\phi$ to make it decay at infinity.

Answer 2

The density that you are looking for is not elementary. To simplify I am taking $k=1.$ You are interested in the characteristic function $\left(\frac{1-iuc}{1-iu}\right)^a$ where $0<c<1$ and $a>0.$ We write $C=(1-c)/c.$ For short let us take $s=iu$ thus dealing with the Laplace transform $$\left(\frac{1-sc}{1-s}\right)^a=c^a\left(1+C\frac{1}{1-s}\right)^a=c^a+c^a\left[\sum_{n=1}^{\infty}\frac{a(a-1)\ldots(a-n+1)}{n!}C^n\frac{1}{(1-s)^n}\right].$$

Now $\int_0^{\infty}e^{sx-x}\frac{x^{n-1}}{(n-1)!}=\frac{1}{(1-s)^n}.$ Therefore consider the function on $R^+$ defined by $$f_a(x)=e^{-x}\left[\sum_{n=1}^{\infty}\frac{a(a-1)\ldots(a-n+1)}{n!}C^n\frac{x^{n-1}}{(n-1)!}\right].$$ As a consequence the desired probability is $$\mu(dx)=c^a\delta_0(dx)+c^af_a(x)1_{R^+}(x)dx.$$ Remark: it is not immediately clear that $f_a$ is positive when $a$ is not an integer. In other terms, it is not clear that $\frac{1-sc}{1-s}$ is the Laplace transform of an infinitely divisible distribution. A short proof of this is using $$\log \frac{1-sc}{1-s}=\log c+\log\left(1+C\frac{1}{1-s}\right)=\log c+\log\left(\sum_{n=1}^{\infty}\frac{(-1)^{n}}{n}C^n\frac{1}{(1-s)^n}\right)=\log c+\int_0^{\infty}e^{-x(1-s)}\frac{1-e^{-Cx}}{x}dx$$ which proves the infinite divisibility with the L'evy measure $e^{-x}\frac{1-e^{-Cx}}{x}1_{R^+}(x)dx.$