I want to reproduce the coefficient estimate of the probit model from the statsmodels function by writing a function that would return the (-loglikelihood) of the probit (standard normal cdf) and the optimize it and return the best iteration. The result should be the same as using the Probit function from statsmodels.
This is the result from the statsmodels, I just want to replicate the estimate of the coefficient and compute the standard errors. Probit function output
data_07["empstat_b"] = data_07["empstat"] == "Employed" # convert empstat to binary
data_07["empstat_b"] = data_07.empstat_b.replace(to_replace=[True, False], value=[1, 0]) # convert to 0,1
Y = data_07["empstat_b"]
X = data_07[['age']]
X = sm.add_constant(X)
model = Probit(Y, X)
probit_model = model.fit()
print(probit_model.summary())