<jupyter_start><jupyter_text>First Method of Simulated Moments (MSM) estimation with estimagicThis tutorial shows how to do a Method of Simulated Moments estimation in estimagic. The Method of Simulated Moments (MSM) is a nonlinear estimation principle that is very useful to fit complicated models to data. The only thing that is needed is a function that simulates model outcomes that you observe in some empirical dataset. In this tutorial we will use a simple linear regression model. This is the same model which we use in the tutorial on maximum likelihood estimation.Throughout the tutorial we only talk about MSM estimation, however, the more general case of indirect inference estimation works exactly the same way.  Steps of MSM estimation- load (simulate) empirical data - define a function to calculate estimation moments on the data - calculate the covariance matrix of the empirical moments (with ``get_moments_cov``)- define a function to simulate moments from the model - estimate the model, calculate standard errors, do sensitivity analysis (with ``estimate_msm``) Example: Estimating the parameters of a regression modelThe model we consider is simple regression model with one variable. The goal is to estimate the slope coefficients and the error variance from a simulated data set.The estimation mechanics are exactly the same for more complicated models. A model is always defined by a function that can take parameters (here: the mean, variance and lower_cutoff and upper_cutoff) and returns a number of simulated moments (mean, variance, soft_min and soft_max of simulated exam points). Model:$$ y = \beta_0 + \beta_1 x + \epsilon, \text{ where } \epsilon \sim N(0, \sigma^2)$$We aim to estimate $\beta_0, \beta_1, \sigma^2$.<jupyter_code>import numpy as np
import pandas as pd

np.random.seed(0)<jupyter_output><empty_output><jupyter_text>Simulate data<jupyter_code>def simulate_data(params, n_draws):
    x = np.random.normal(0, 1, size=n_draws)
    e = np.random.normal(0, params.loc["sd", "value"], size=n_draws)
    y = params.loc["intercept", "value"] + params.loc["slope", "value"] * x + e
    return pd.DataFrame({"y": y, "x": x})
true_params = pd.DataFrame(
    data=[[2, -np.inf], [-1, -np.inf], [1, 1e-10]],
    columns=["value", "lower_bound"],
    index=["intercept", "slope", "sd"],
)

data = simulate_data(true_params, n_draws=100)<jupyter_output><empty_output><jupyter_text>Calculate Moments<jupyter_code>def calculate_moments(sample):
    moments = {
        "y_mean": sample["y"].mean(),
        "x_mean": sample["x"].mean(),
        "yx_mean": (sample["y"] * sample["x"]).mean(),
        "y_sqrd_mean": (sample["y"] ** 2).mean(),
        "x_sqrd_mean": (sample["x"] ** 2).mean(),
    }
    return pd.Series(moments)
empirical_moments = calculate_moments(data)
empirical_moments<jupyter_output><empty_output><jupyter_text>Calculate the covariance matrix of empirical momentsThe covariance matrix of the empirical moments (``moments_cov``) is needed for three things:1. to calculate the weighting matrix2. to calculate standard errors3. to calculate sensitivity measuresWe will calculate ``moments_cov`` via a bootstrap. Depending on your problem there can be other ways to do it.<jupyter_code>from estimagic import get_moments_cov
moments_cov = get_moments_cov(
    data, calculate_moments, bootstrap_kwargs={"n_draws": 5_000, "seed": 0}
)

moments_cov<jupyter_output><empty_output><jupyter_text>``get_moments_cov`` mainly just calls estimagic's bootstrap function. See our [bootstrap_tutorial](../how_to_guides/inference/how_to_do_bootstrap_inference.ipynb) for background information.   Define a function to calculate simulated momentsIn a real application, this is the step that takes most of the time. However, in our very simple example, all the work is already done by numpy.<jupyter_code>def simulate_moments(params, n_draws=10_000, seed=0):
    np.random.seed(seed)
    sim_data = simulate_data(params, n_draws)
    sim_moments = calculate_moments(sim_data)
    return sim_moments
simulate_moments(true_params)<jupyter_output><empty_output><jupyter_text>Estimate the model parametersEstimating a model consists of the following steps:- Building a criterion function that measures a distance between simulated and empirical moments- Minimizing this criterion function- Calculating the Jacobian of the model- Calculating standard errors, confidence intervals and p values- Calculating sensitivity measuresThis can all be done in one go with the ``estimate_msm`` function. This function has good default values, so you only need a minimum number of inputs. However, you can configure almost every aspect of the workflow via optional arguments. If you need even more control, you can call the low level functions ``estimate_msm`` is built on directly.<jupyter_code>from estimagic import estimate_msm
start_params = true_params.assign(value=[100, 100, 100])

res = estimate_msm(
    simulate_moments,
    empirical_moments,
    moments_cov,
    start_params,
    optimize_options={"algorithm": "scipy_lbfgsb"},
)
res["summary"]<jupyter_output><empty_output><jupyter_text>What's in the result?The result is a dictionary with the following entries:<jupyter_code>res.keys()<jupyter_output><empty_output><jupyter_text>- summary: A DataFrame with estimated parameters, standard errors, p_values and confidence intervals. This is ideal for use with our `estimation_table` function. - cov: A DataFrame with the full covariance matrix of the estimated parameters. You can use this for hypothesis tests. - optimize_res: A dictionary with the complete output of the numerical optimization (if one was performed)- numdiff_info: A dictionary with the complete output of the numerical differentiation (if one was performed)- jacobian: A DataFrame with the jacobian of `simulate_moments` with respect to the free parameters. - sensitivity: A dictionary of DataFrames with the six sensitivity measures from Honore, Jorgensen and DePaula (2020)  How to visualize sensitivity measuresFor more background on the sensitivity measures and their interpretation, check out the [how to guide](../how_to_guides/miscellaneous/how_to_visualize_and_interpret_sensitivity_measures.ipynb) on sensitivity measures. Here we are just showing you how to plot them:<jupyter_code>from estimagic.visualization.lollipop_plot import lollipop_plot

sensitivity_data = res["sensitivity"]["sensitivity_to_bias"].abs().T

fig = lollipop_plot(sensitivity_data)

fig = fig.update_layout(height=500, width=900)
fig.show(renderer="png")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08285105086810844, 'fraction_numerical': 0.00837648492232714, 'mean_word_length': 5.6602434077079105, 'pattern_counts': {'":': 10, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '148904', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1938, 'n_tokens_neox': 1822, 'n_words': 718, 'path': 'docs/source/getting_started/first_msm_estimation_with_estimagic.ipynb', 'repo_name': 'janosg/estimagic', 'size': '6566'}