Numerical Inverse Z-Transform - Abate and Whitt

Question

I'm trying to implement an inverse Z-Transform using the Fourier series based technique by Abate and Whitt:

The Fourier-series method for inverting transforms of probability distributions.

Numerical inversion of probability generating functions.

We approximate the integral numerically by $$f(n) \approx \frac{1}{2n\rho^n}\bigg[ \tilde{f}(\rho) + 2 \sum^{n-1}_{j = 1} (-1)^j\text{Re}\tilde{f}(\rho e^{\frac{j\pi i}{n}}) + (-1)^n\tilde{f}(-\rho) \bigg]$$

and setting $\rho = 10^{-\lambda/2n}$ gives an accuracy of $10^{-\lambda}$.

import cmath
import numpy as np

def abate_whitt(ftilde, n):
    lmbda = 1
    rho = 10 ** (-lmbda / (2 * n))

    e = cmath.e
    pi = cmath.pi

    summation = 0
    for k in range(1, n):
        summation += ((-1)**k) * (np.real(ftilde(rho * e**((1j * pi * k) / n))))
    summation *= 2
    part1 = ftilde(rho)
    part2 = ((-1)**n) * ftilde(-rho)
    scale_factor = 1 / (2 * n * (rho**n))

    result = scale_factor * (part1 + summation + part2)
    return result

I then test it for pairs of functions and started off with $\tilde{f}(z) = \frac{z}{z-1}$ and $f(t) = 1$.

def ftilde(z):
   return z / (z-1)

Now from my understanding, given any values of n, the result should be 1. However, given any n values, I get around -1.111 (3 d.p.). Would anyone be able to help me understand why?

Full Code + Result

Answer 1

Can you please tell me where you found the above formula?
If you have poles on or outside the unit circle you may use this formula:

$$f(\tilde{f},T,\rho,n)\approx \frac{\rho ^T}{n}\cdot \sum _{k=0}^{n-1} \Re\left(\tilde{f}\left(\rho\cdot \exp \left(\frac{i ((2 \pi ) k)}{n}\right)\right)\cdot \exp \left(\frac{i ((2 \pi ) k T)}{n}\right) \right)$$

IMPORTANT: $\tilde{f}(z)$ all poles are located inside the circle with radius $\rho$ and has no poles at $\infty$.

import matplotlib.pyplot as plt
import numpy as np

# Numerical approximation of inverse Z-Transform.
def abate_whitt(ftilde, T, rho=1, n=50):
  sum = 0.0
  for k in range(0, n):
    exp1 = np.exp(1j * 2*np.pi * k / n)
    exp2 = exp1**T
    sum += ftilde(rho*exp1) * exp2
  return (rho**T / n) * sum.real

# Z-functions for test.
def ftilde1(z): return z / (z - 1)
def ftilde2(z): return z / ((z + 1/2) * (z + 1/3))

# Exact solution
def f1(T): return 1
def f2(T): return 6*(-(-1/2)**T + (-1/3)**T)

T    = 1
val1 = abate_whitt(ftilde1, T, 2)
val2 = f1(T)
print("value=", val1, "exact=", val2, "error=", val1-val2)

value= 1.0000000000000007 exact= 1 error= 6.661338147750939e-16

Preferred inverse Z-Transform with IFFT

n   = 80
k   = np.arange(n)
phi = 2 * np.pi * k/n
F   = ftilde2(np.exp(1j * phi))
fapp= np.fft.ifft(F).real

plt.figure()
plt.plot(k, fapp-f2(k))
plt.title('Approximation error')