Pandas is calculating the UNBIASED estimator of the excess Kurtosis. Kurtosis is the normalized 4th central moment. To find the unbiased estimators of the cumulants you need the k-statistics.

So the unbiased estimator of kurtosis is (k4/k2**2)

To illustrate this:

import pandas as pd
import numpy as np

np.random.seed(11234)
test_series = pd.Series(np.random.randn(5000))
test_series.kurtosis()
#-0.0411811269445872

Now we can calculate this explicitly using the k-statistics:

n = len(test_series)
S1 = test_series.pow(1).sum()
S2 = test_series.pow(2).sum()
S3 = test_series.pow(3).sum()
S4 = test_series.pow(4).sum()

# Eq (7) and (5) from the k-statistics link
k4 = (-6*S1**4 + 12*n*S1**2*S2 - 3*n*(n-1)*S2**2 -4*n*(n+1)*S1*S3 + n**2*(n+1)*S4)/(n*(n-1)*(n-2)*(n-3))
k2 = (n*S2-S1**2)/(n*(n-1))
# k2 is the same as the N-1 variance: test_series.std(ddof=1)**2

k4/k2**2
#-0.04118112694458816

If you want better agreement to more decimal places, you'll need to be careful with the sums as they get rather large. But they're identical to 12 places.

Pandas calculate UNBIASED estimator of the population kurtosis. Look at the Wikipedia for formulas: https://www.wikiwand.com/en/Kurtosis

Calculate kurtosis from scratch

import numpy as np
import pandas as pd
import scipy

x = np.array([0, 3, 4, 1, 2, 3, 0, 2, 1, 3, 2, 0,
              2, 2, 3, 2, 5, 2, 3, 999])
xbar = np.mean(x)
n = x.size
k2 = x.var(ddof=1) # default numpy is biased, ddof = 0
sum_term = ((x-xbar)**4).sum()
factor = (n+1) * n / (n-1) / (n-2) / (n-3)
second = - 3 * (n-1) * (n-1) / (n-2) / (n-3)

first = factor * sum_term / k2 / k2

G2 = first + second
G2 # 19.998428728659768

Calculate kurtosis using numpy/scipy

scipy.stats.kurtosis(x,bias=False) # 19.998428728659757

Calculate kurtosis using pandas

pd.DataFrame(x).kurtosis() # 19.998429

Similarly, you can also calculate skewness.