Download the notebook here

Method of Simulated Moments (MSM)

This 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 model

The 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\).

import numpy as np
import pandas as pd

np.random.seed(0)

Simulate data

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)

Calculate Moments

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

y_mean         2.022205
x_mean         0.059808
yx_mean       -0.778369
y_sqrd_mean    5.942648
x_sqrd_mean    1.019404
dtype: float64

Calculate the covariance matrix of empirical moments

The covariance matrix of the empirical moments (moments_cov) is needed for three things: 1. to calculate the weighting matrix 2. to calculate standard errors 3. to calculate sensitivity measures

We will calculate moments_cov via a bootstrap. Depending on your problem there can be other ways to do it.

from estimagic import get_moments_cov

moments_cov = get_moments_cov(
    data, calculate_moments, bootstrap_kwargs={"n_draws": 5_000, "seed": 0}
)

moments_cov

	y_mean	x_mean	yx_mean	y_sqrd_mean	x_sqrd_mean
y_mean	0.019384	-0.009285	-0.014657	0.076048	-0.002489
x_mean	-0.009285	0.010238	0.018576	-0.035035	0.001544
yx_mean	-0.014657	0.018576	0.054319	-0.081187	-0.010329
y_sqrd_mean	0.076048	-0.035035	-0.081187	0.353938	0.000089
x_sqrd_mean	-0.002489	0.001544	-0.010329	0.000089	0.015506

get_moments_cov mainly just calls estimagic’s bootstrap function. See our bootstrap_tutorial for background information.

Define a function to calculate simulated moments

In 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.

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)

y_mean         2.029422
x_mean        -0.018434
yx_mean       -1.020287
y_sqrd_mean    6.095197
x_sqrd_mean    0.975608
dtype: float64

Estimate the model parameters

Estimating 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 measures

This 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.

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="scipy_lbfgsb",
)

res.summary()

	value	standard_error	p_value	ci_lower	ci_upper	stars
intercept	1.995794	0.138153	2.645631e-47	1.725019	2.266568	***
slope	-0.750961	0.238691	1.654246e-03	-1.218786	-0.283136	***
sd	1.143843	0.149262	1.811586e-14	0.851296	1.436391	***

What’s in the result?

MomentsResult objects provide attributes and methods to calculate standard errors, confidence intervals and p_values with multiple methods. You can even calculate cluster robust standard errors.

A few examples are:

res.params

	value	lower_bound
intercept	1.995794	-inf
slope	-0.750961	-inf
sd	1.143843	1.000000e-10

res.cov(method="robust")

	intercept	slope	sd
intercept	0.019086	-0.013996	-0.010416
slope	-0.013996	0.056973	0.027351
sd	-0.010416	0.027351	0.022279

res.se()

	value	lower_bound
intercept	0.138153	-inf
slope	0.238691	-inf
sd	0.149262	1.000000e-10

How to visualize sensitivity measures

For more background on the sensitivity measures and their interpretation, check out the how to guide on sensitivity measures.

Here we are just showing you how to plot them:

from estimagic.visualization.lollipop_plot import lollipop_plot

sensitivity_data = res.sensitivity(kind="bias").abs().T

fig = lollipop_plot(sensitivity_data)

fig = fig.update_layout(height=500, width=900)
fig.show(renderer="png")

previous

Likelihood estimation

next

Numerical differentiation