Why doing bootstrap with `sklearn` and doing bootstrap with `statsmodels` yields different coefficients?

I have a dataset: link to download here.

I'm performing LogisticRegression using both sklearn and statsmodels.

When doing bootstrap to estimate coef and std_err, the coef yielded from the bootstrap_sklearn is different from the coef that I get doing non_bootstrap_sklearn, non_bootstrap_statsmodels and bootstrap_statsmodels.

First, the 'real' coefficients using statsmodels:

endog = (Default["default"] == "Yes").astype(int)
exog = sm.add_constant(Default[["income", "balance"]])
mod = sm.Logit(endog, exog)
res = mod.fit()
res.summary()

Output:

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        -11.5405      0.435    -26.544      0.000     -12.393     -10.688
income      2.081e-05   4.99e-06      4.174      0.000     1.1e-05    3.06e-05
balance        0.0056    0.00022     24.835      0.000       0.005       0.006
==============================================================================

We have income_coef = 2e-05, income_std_err = 5e-06, balance_coef = 0.0056 and balance_std_err = 0.00022.

I then perform bootstrap using sklearn.LogisticRegression:

np.random.seed(17)
num_estimates = 1000
boot_estimates = np.empty((num_estimates, 2))

for i in range(num_estimates):
    clf = LogisticRegression(penalty = "none", solver = "lbfgs")
    sample = resample(Default)
    coefs = boot_fn(sample, clf)
    boot_estimates[i, 0] = coefs[0, 0]
    boot_estimates[i, 1] = coefs[0, 1]
    print(f"{i+1} out of {num_estimates}", end="\r")

And bootstrap with statsmodels.Logit:

np.random.seed(17)
num_estimates = 1000
boot_estimates = np.empty((num_estimates, 2))
for i in range(num_estimates):
    sample = resample(Default)
    endog = (sample["default"] == "Yes").astype(int)
    exog = sm.add_constant(sample[["income", "balance"]])
    mod = sm.Logit(endog, exog)
    res = mod.fit(disp = False)
    boot_estimates[i, 0] = res.params["income"]
    boot_estimates[i, 1] = res.params["balance"]
    print(f'{i+1} of {num_estimates}', end='\r')

I then add everything to a df for easy comparison (using .mean() and .sem() to get coef and std_err), here's the df:

                           income   balance
statsmodels_coef     2.080898e-05  0.005647
bs_sk_coef          -6.101152e-05  0.002728
statsmodels_std_err  4.985245e-06  0.000227
bs_sk_std_err        2.322953e-06  0.000082
bs_sm_coef           2.072321e-05  0.005655
bs_sm_std_err        1.510824e-07  0.000007

-------------EDIT--------------

Here's the definition of boot_fn:

def boot_fn(data, clf):
    clf.fit(data[["income", "balance"]], data["default"])
    return clf.coef_

Furthermore, coef_ with sklearn:

X = Default[["balance", "income"]]
y = Default["default"]
clf = LogisticRegression(penalty = "none", solver = "lbfgs")
clf.fit(X, y)
clf.coef_

Output:

array([[5.64710291e-03, 2.08089921e-05]])

To conclude, the 'real' statsmodels_coef and sklearn_coef are the same, together with the bootstrap_statsmodels_coef but the bootstrap_sklearn_coef yields a different result.

Why is that?