<jupyter_start><jupyter_text>First maximum likelihood estimation with estimatigThis notebook shows how to do a simple maximum likelihood (ml) estimation with estimagic. As illustrating example we implement an ordered logit model from scratch. 1. Load and process the data2. Set up a likelihood function3. Maximize the likelihood function4. Calculate standard errors, P-values and confidence intervalsThe user only needs to do step 1 and 2. The rest is done by `estimate_ml`. To be very clear: Estimagic is not a package to estimate logit models or other models that are implemented in Stata, statsmodels or anywhere else. Its purpose is to estimate parameters with custom likelihood or method of simulated moments functions. We just use an ordered logit model as an example of a very simple likelihood function. The example we will use to test our model is taken from the [Stata Documentation](https://stats.idre.ucla.edu/stata/dae/ordered-logistic-regression/).<jupyter_code>import numpy as np
import pandas as pd
from patsy import dmatrices
from scipy import stats

from estimagic import estimate_ml<jupyter_output><empty_output><jupyter_text>Process the user inputWe choose an R-style formula as a convenient way of specifying the ordered logit model and use `patsy` to construct matrices from the dataset. We will need four inputs:1. A DataFrame with start parameters for the optimization.2. An array with the dependent variable.3. A 2d array with explanatory variables.4. Constraints for the optimization that keep the cutoffs increasing.We construct all of those inputs using the `ordered_logit_processing` function. You could also do those steps in a simple script.<jupyter_code>def ordered_logit_processing(formula, data):
    """Process user input for an ordered logit model."""
    # extract data arrays
    y, x = dmatrices(formula + " - 1", data, return_type="dataframe")
    y = y[y.columns[0]]

    # extract dimensions
    num_choices = len(y.unique())
    beta_names = list(x.columns)
    num_betas = len(beta_names)
    num_cutoffs = num_choices - 1

    # set-up index for params_df
    names = beta_names + list(range(num_cutoffs))
    categories = ["beta"] * num_betas + ["cutoff"] * num_cutoffs
    index = pd.MultiIndex.from_tuples(zip(categories, names), names=["type", "name"])

    # make params_df
    np.random.seed(5471)
    start_params = pd.DataFrame(index=index)
    start_params["value"] = np.hstack(
        [
            np.random.uniform(low=-0.5, high=0.5, size=len(x.columns)),
            np.arange(num_cutoffs) * 2,
        ]
    )
    start_params["group"] = start_params.index.get_level_values("type")

    # make constraints
    constr = [{"loc": "cutoff", "type": "increasing"}]

    # turn pandas objects into numpy arrays
    y_arr = y.to_numpy().astype(int)
    x_arr = x.to_numpy()

    return start_params, y_arr, x_arr, constr<jupyter_output><empty_output><jupyter_text>Defining the `loglike` function<jupyter_code>def ordered_logit_loglike(params, y, x):
    """Likelihood function of an orderd logit model."""
    # parse the parameter vector into its quantities
    beta = params.loc["beta", "value"].to_numpy()
    cutoffs = params.loc["cutoff", "value"].to_numpy()

    # calculate deterministic part of utilities
    xb = x.dot(beta)

    # evaluate likelihood
    upper_cutoffs = np.hstack([cutoffs, np.inf])[y]
    lower_cutoffs = np.hstack([-np.inf, cutoffs])[y]
    upper_cdf = stats.logistic.cdf(upper_cutoffs - xb)
    lower_cdf = stats.logistic.cdf(lower_cutoffs - xb)

    contributions = np.log(upper_cdf - lower_cdf)

    res = {"contributions": contributions, "value": contributions.sum()}

    return res<jupyter_output><empty_output><jupyter_text>A few remarks are in order:1. There are numerically better ways to calculate the likelihood, we chose this implementation for brevity and readability. 2. The loglike function takes params and other arguments. You are completely flexible in the number and names of the other arguments as long as the first argument is params. 3. The loglike function returns a dictionary with the entries "contributions" and "value". The "contributions" are the log likelihood evaluations of each individual in the dataset. The "value" are their sum. The "value" entry could be omitted, the "contributions" entry is mandatory.   Estimating the model<jupyter_code>data = pd.read_pickle("ologit.pickle")
formula = "apply ~ pared + public + gpa"
start_params, y, x, constraints = ordered_logit_processing(formula, data)

res = estimate_ml(
    loglike=ordered_logit_loglike,
    params=start_params,
    optimize_options={"algorithm": "scipy_lbfgsb"},
    constraints=constraints,
    loglike_kwargs={"y": y, "x": x},
)
res["summary_jacobian"].round(3)<jupyter_output><empty_output><jupyter_text>What's in the results?The result of `estimate_ml` is a dictionary with the following entries:<jupyter_code>res.keys()<jupyter_output><empty_output><jupyter_text>Importantly, we might have several summaries and several cov entries. This is because we always all possible types of standard errors, so you can compare them easily (for example using our estimation table functions). In the current example we only have jacobian based standard errors because we did not provide a closed form hessian function and numerical hessians are not yet implemented. With a closed form hessian we would also get hessian based and robust standard errors. Even cluster and strata robust standard errors are possible if you provide the relevant information. If numerical optimizations or derivative calculations were performed, the full output of those steps is also part of the results dictionary.   Compare to STATA's results<jupyter_code>stata_results = pd.read_csv("stata_ologit_results.csv")
stata_results.round(3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.082116474832503, 'fraction_numerical': 0.004810170073870469, 'mean_word_length': 4.685546875, 'pattern_counts': {'":': 7, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '296047', 'license': 'MIT', 'n_tokens_mistral': 1699, 'n_tokens_neox': 1644, 'n_words': 630, 'path': 'docs/source/getting_started/first_likelihood_estimation_with_estimagic.ipynb', 'repo_name': 'OpenSourceEconomics/estimagic', 'size': '5821'}