Estimating truncation error in hypergeometric series

Question

If $\sum_{k=0}^\infty s_k$ is a hypergeometric series, $g(n)=\frac{s_{n}}{s_{n-1}}$ is a simple expression which can be written $g(n)=\sum_{k=0}^\infty \frac{b_k}{n^k}$.

I want to estimate the truncation error $e(n)=\frac{1}{s_n}\cdot\sum_{k=n}^\infty s_k=\sum_{k=0}^\infty \frac{a_k}{n^k}$.

From $\sum_{k=n-1}^\infty s_k = s_{n-1} + \sum_{k=n}^\infty s_k$ we obtain $$e(n-1)=1+e(n)\cdot g(n).$$ This equation can be used to calculate the coefficients $a_k$ recursively: The constant terms yield $a_0 = \frac{1}{1-b_0}$, and for $N\geq 1$ we have $$a_N = \frac{b_N}{(1-b_0)^2}-\sum_{k=1}^{N-1} a_k\cdot \frac{\binom{N-1}{k-1}-b_{N-k}}{1-b_0}$$

My Question:

Writing $e(n)=\sum_{k=0}^3 \frac{a_k}{n^k} + R_4(n)$, how can I get bounds on $R_4(n)$?

Numerical Example: For the Chudnovsky Formula, the following python-code yields the values of $b_k$ and $a_k$. The results show that $a_k$ is mostly alternating, only $a_0$ and $a_1$ have the same sign, and $a_{47}$ and $a_{48}$. From $a_{48}$ on, the absolute values of $a_n$ seem to be growing exponentially.

MAXORDER=70

from fractions import Fraction as frac
# def frac(a,b): return float(a)/float(b)

# Pascal's Triangle from https://stackoverflow.com/a/24093559/11225178
def pascal(nrows):
    results = []
    for _ in range(nrows+1): 
        row = [1]
        if results:
            lastrow = results[-1]
            row.extend([sum(pair) for pair in zip(lastrow, lastrow[1:])])
            row.append(1)
        results.append(row)
    return results

binom=pascal(MAXORDER)

# calculate b_k=B[k] for the Chudnovsky Algorithm:
y = frac(1,1)-frac(13591409,545140134)
z = frac(-1,53360**3)
B = [z*frac(1,1), z*frac(-1,2), z*(y-frac(31,36)), z*(y**2-y*frac(3,2)+frac(41,72)), z*(y**3-frac(3,2)*y**2+frac(23,36)*y-frac(5,72))]
while (len(B)<MAXORDER+1): B.append(B[-1]*y)

# calculate a_k=A[k] for the Chudnovsky Algorithm:
A = [1/(1-z)]
while (len(A)<len(B)):
  N = len(A)
  k = 1
  sum = B[N]/(1-z)
  while (k<N) :
    sum -= A[k]*(binom[N-1][k-1]-B[N-k])
    k = k+1
  A.append(sum/(1-z))

print("values of k, b_k, a_k below:")
k = 0
while k<len(A):
    print(k, float(B[k]), float(A[k]))
    k=k+1

EDIT: Formula for Chudnovsky-Example In case of the Chudnovsky-Series, we have $$s_n = \frac{(-1)^n\left(6n\right)!}{\left(3n\right)!\left(n!\right)^3}\,\frac{S + n}{640320^{3n}}$$ with $S=\frac{13591409}{545140134}$ and thus $$g(n)=\frac{s_{n}}{s_{n-1}} =-53360^{-3}\cdot\left(1-\frac{3}{2n}+\frac{23}{36n^2}-\frac{5}{72n^3}\right)\cdot \left(1+\frac{1}{S+n-1}\right).$$

EDIT: Divergent series The following graph shows that $-\frac{a_k}{a_{k-1}}\approx \frac{k}{33}$ and thus $|a_k|\approx |a_1|\cdot\frac{k!}{33^k}$ which proves $e(n)$ to be a divergent series.

But still, $\sum_{k=0}^3 \frac{a_k}{n^k}$ is quite a good approximation for $e(n)$.