Method of Simulated Moments Estimation (MSM) — EKW lectures documentation (2025)

Introduction

The method of simulated moments approach to estimating model parameters is to minimize a certain distance between observed moments and simulated moments with respect to the parameters that generate the simulated model.

Denote \(X\) our observed data and \(m(X)\) the vector of observed moments. To construct the criterion function, we use the parameter vector \(\theta\) to simulate data from the model \(\hat{X}\). We can then calculate the simulated moments \(m(\hat{X}| \theta)\).

The criterion function is then given by

\begin{equation}\psi(\theta) = (m(X) - m(\hat{X}| \theta))'\Omega(m(X) - m(\hat{X}| \theta))\end{equation}

where the difference between observed and simulated moments \((m(X) - m(\hat{X}| \theta))\) constitutes a vector of the dimension \(1\times M\) with \(1,..,M\) denoting the number of moments. The \(M\times M\) weighting matrix is given by \(\Omega\).

The MSM estimator is defined as the solution to

\begin{equation}\hat{\theta}(\Omega) = \underset{\theta}{\operatorname{argmin}} \psi(\theta).\end{equation}

The criterion function is thus a strictly positive scalar and the estimator depends on the choice of moments \(m\) and the weighting matrix \(\Omega\). The weighting matrix applies some scaling for discrepancies between the observed and simulated moments. If we use the identity matrix, each moment is given equal weight and the criterion function reduces to the sum of squared moment deviations.

Aside from the choice of moments and weighting matrix, some other important choices that influence the the estimation are the simulator itself and the algorithm and specifications for the optimization procedure. Many explanations for simulated method of moments estimation also feature the number of simulations as a factor that is to be determined for estimation. We can ignore this factor for now since we are working with a large simulated dataset.

In the following we will set up the data, moments and weighting matrix needed to construct the criterion function and subsequently try to estimating the parameters of the model using this MSM setup.

Observed Data

We generate our model and data using respy and a simple Robinson Crusoe model. In this model, the agent, Robinson Crusoe, in each time period decides between two choice options: working (i.e. going fishing) or spending time in the hammock.

We can use respy to simulate the data for this exercise.

[2]:

params_true, options = load_model_specs()

[3]:

# Generate observed data from model.simulate = rp.get_simulate_func(params_true, options)data_obs = simulate(params_true)

Let’s take a look at the model specifications.

[4]:

params_true

[4]:

[5]:

options

[5]:

{'estimation_draws': 100, 'estimation_seed': 100, 'estimation_tau': 0.001, 'interpolation_points': -1, 'n_periods': 5, 'simulation_agents': 1000, 'simulation_seed': 132, 'solution_draws': 100, 'solution_seed': 456, 'covariates': {'constant': '1'}}

[6]:

data_obs.head(10)

[6]:

Data Moments

For the setup of the estimation we first have to choose a set of moments that we will use to match the observed data and the simulated model. For this model we include two sets of moments:

1. The first set are Robinson’s choice frequencies (choice frequencies here refers to the share of agents that have chosen a specific option) for each period.

2. The second set are moments that characterize the wage distribution for each period, i.e. the mean of the wage of all agents that have chosen fishing in a given period and the standard deviation of the wages.

We need a functions that compute the set of moments on the observed and simulated data.

[7]:

def calc_choice_frequencies(df): """Calculate choice frequencies""" return df.groupby("Period").Choice.value_counts(normalize=True).unstack()

[8]:

def calc_wage_distribution(df): """Calculate wage distribution.""" return df.groupby(["Period"])["Wage"].describe()[["mean", "std"]]

[9]:

# Save functions in a dictionary to pass to the criterion function later on.calc_moments = { "Choice Frequencies": calc_choice_frequencies, "Wage Distribution": calc_wage_distribution,}

Respy’s MSM interface additionally requires us to specify, how missings in the computed moments should be handled. Missing moments for instance can arise, when certain choice options aren’t picked at all for some periods. We just specify a function that replaces missings with 0.

[10]:

def replace_nans(df): """Replace missing values in data.""" return df.fillna(0)

Now we are ready to calculate the moments.

[11]:

moments_obs = { "Choice Frequencies": replace_nans(calc_moments["Choice Frequencies"](data_obs)), "Wage Distribution": replace_nans(calc_moments["Wage Distribution"](data_obs)),}

[12]:

print("Choice Frequencies")print(moments_obs["Choice Frequencies"])print("\n Wage Distribution")print(moments_obs["Wage Distribution"])

Choice Frequencies fishing hammockPeriod0 0.698 0.3021 0.698 0.3022 0.698 0.3023 0.698 0.3024 0.698 0.302 Wage Distribution mean stdPeriod0 1.003189 0.0086401 1.072815 0.0107372 1.150106 0.0123163 1.234041 0.0124064 1.322711 0.013473

Weighting Matrix

Next we specify a weighting matrix. It needs to be a square matrix with the same number of diagonal elements as there are moments. One option would be to use the identity matrix, but we use a weighting matrix that adjusts for the variance of each moment to improve the estimation process. The variances for the moments are constructed using a bootstrapping procedure.

[13]:

def get_weighting_matrix(data, calc_moments, num_boots, num_agents_msm): """ Compute weighting matrix for estimation with MSM.""" # Seed for reproducibility. np.random.seed(123) index_base = data.index.get_level_values("Identifier").unique() # Create bootstrapped moments. moments_sample = list() for _ in range(num_boots): ids_boot = np.random.choice(index_base, num_agents_msm, replace=False) moments_boot = [calc_moments[key](data.loc[ids_boot, :]) for key in calc_moments.keys()] moments_boot = rp.get_flat_moments(moments_boot) moments_sample.append(moments_boot) # Compute variance for each moment and construct diagonal weighting matrix. moments_var = np.array(moments_sample).var(axis=0) weighting_matrix = np.diag(moments_var ** (-1)) return np.nan_to_num(weighting_matrix)

[14]:

W = get_weighting_matrix(data_obs, calc_moments, 300, 500)

[15]:

pd.DataFrame(W)

[15]:

Criterion Function

We have collected the observed data for our model, chosen the set of moments we want to use for estimation and defined a weighting matrix based on these moments. We can now set up the criterion function to use for estimation.

As already discussed above, the criterion function is given by the weighted square product of the difference between observed moments \(m(X)\) and simulated moments \(m(\hat{X}| \theta)\). Trivially, if we have that \(m(X) = m(\hat{X}| \theta)\), the criterion function returns a value of 0. Thus, the closer \(\theta\) is to the real parameter vector, the smaller should be the value for the criterion function.

[16]:

criterion_msm = rp.get_moment_errors_func( params_true, options, calc_moments, replace_nans, moments_obs, W)

Criterion function at the true parameter vector:

[17]:

fval = criterion_msm(params_true)fval

[17]:

0.0

We can plot the criterion function to examine its behavior around the minimum in more detail. The plots below show the criterion function at varying values of all parameters in the the paramter vector.

[18]:

plot_criterion_params(params_true, list(params_true.index), criterion_msm, 0.05)

This depiction for most parameters conceals the fact that the criterion function is not a smooth function of our parameter values. We can reveal this property if we ‘zoom in’ on the function far enough. The plots below show the criterion function for varying values of delta around the true minimum value of 0.95. We can see that the function exhibits small plateaus and is thus not completely smooth.

[19]:

for radius in [0.05, 0.01, 0.001, 0.0001]: plot_criterion_params(params_true, list(params_true.index)[0:1], criterion_msm, radius)

Estimation Exercise

In the following we will conduct a simulation exercise to estimate the parameter vector using our criterion function and weighting matrix. We will begin by simulating data using the new parameter vector and examine how the simulated moments differ from the observed ones. We will then use an optimizer to minimize the criterion function in order to retrieve the true parameter vector. Additionally, we will explore how the criterion function behaves if we change the simulation seed, move away fromthe true values and test the optimization for a derivative-based optimization algorithm.

[20]:

params_cand = params_true.copy()params_cand.loc["delta", "value"] = 0.93

[21]:

params = params_cand.copy()simulate = rp.get_simulate_func(params, options)df_sim = simulate(params)

[22]:

moments_sim = { "Choice Frequencies": replace_nans(calc_moments["Choice Frequencies"](df_sim)), "Wage Distribution": replace_nans(calc_moments["Wage Distribution"](df_sim)),}

[23]:

plot_moments_choices(moments_obs, moments_sim)

[24]:

plot_moments_wage(moments_obs, moments_sim)

[25]:

fval = criterion_msm(params_cand)fval

[25]:

7838.131228444257

[26]:

from estimagic import minimize # noqa: E402

[27]:

rslt = minimize( criterion=criterion_msm, params=params_true, algorithm="nag_pybobyqa",)rslt["solution_criterion"]

[27]:

0.0

[28]:

params_cand["lower_bound"] = [0.69, 0.066, -0.2, 1, 0, 0, 0]params_cand["upper_bound"] = [1, 0.078, -0.09, 1.1, 0.5, 0.5, 0.5]params_cand

[28]:

[29]:

# Define base constraints to use for the rest of the notebook.constr_base = [ {"loc": "shocks_sdcorr", "type": "sdcorr"}, {"loc": "delta", "type": "fixed"}, {"loc": "wage_fishing", "type": "fixed"}, {"loc": "nonpec_fishing", "type": "fixed"}, {"loc": "nonpec_hammock", "type": "fixed"}, {"loc": "shocks_sdcorr", "type": "fixed"},]

[30]:

# Remove constraint for delta.constr = constr_base.copy()constr.remove({"loc": "delta", "type": "fixed"})

[31]:

rslt = minimize( criterion=criterion_msm, params=params_cand, algorithm="nag_pybobyqa", constraints=constr,)rslt["solution_criterion"]

[31]:

0.5248952065747418

[32]:

options_new_seed = options.copy()options_new_seed["simulation_seed"] = 400

[33]:

criterion_msm = rp.get_moment_errors_func( params_true, options_new_seed, calc_moments, replace_nans, moments_obs, W)criterion_msm(params_true)

[33]:

38.97135021944551

[34]:

rslt_new_seed = minimize( criterion=criterion_msm, params=params_true, algorithm="nag_pybobyqa", constraints=constr,)rslt_new_seed["solution_criterion"]

[34]:

17.369549302805428

[35]:

rslt_new_seed = minimize( criterion=criterion_msm, params=params_cand, algorithm="nag_pybobyqa", constraints=constr,)rslt_new_seed["solution_criterion"]

[35]:

17.369549302805428

[36]:

# List of seeds, includes true simulation seed of 132.seeds = list(range(120, 140, 1))

[37]:

# Define other inputs for criterion functioncriterion_kwargs = dict()criterion_kwargs["params"] = params_true.copy()criterion_kwargs["options"] = optionscriterion_kwargs["calc_moments"] = calc_momentscriterion_kwargs["replace_nans"] = replace_nanscriterion_kwargs["moments_obs"] = moments_obscriterion_kwargs["weighting_matrix"] = W

[38]:

plot_chatter(seeds, criterion_kwargs)

[39]:

# Changing the number of agents.num_agents = [500, 1000, 5000, 10000, 50000]

[40]:

plot_chatter_numagents_sim(seeds, num_agents, criterion_kwargs)

[41]:

plot_chatter_numagents_both(seeds, num_agents, calc_moments, replace_nans, criterion_kwargs)

[42]:

params_cand.loc["wage_fishing", "value"] = 0.072params_cand

[42]:

[43]:

criterion_msm = rp.get_moment_errors_func( params_true, options_new_seed, calc_moments, replace_nans, moments_obs, W)rslt_wrong_fix = minimize( criterion=criterion_msm, params=params_cand, algorithm="nag_pybobyqa", constraints=constr,)rslt_wrong_fix["solution_criterion"]

[43]:

813.7988993578208

[44]:

params_cand.loc[:, "value"] = rslt_wrong_fix["solution_params"][["value"]]params_cand

[44]:

[45]:

# Adjust constraints to free up both delta and wage_fishing.constr_u = constr_base.copy()constr_u.remove({"loc": "delta", "type": "fixed"})constr_u.remove({"loc": "wage_fishing", "type": "fixed"})

[62]:

rslt_unfix = minimize( criterion=criterion_msm, params=params_cand[["value"]], algorithm="nag_pybobyqa", constraints=constr_u,)rslt_unfix["solution_criterion"]

[62]:

710.9165575145821

[63]:

deviation = params_true["value"] - rslt_unfix["solution_params"]["value"]deviation

[63]:

category namedelta delta 0.026703wage_fishing exp_fishing -0.001791nonpec_fishing constant 0.000000nonpec_hammock constant 0.000000shocks_sdcorr sd_fishing 0.000000 sd_hammock 0.000000 corr_hammock_fishing 0.000000Name: value, dtype: float64

[64]:

# Define candidate parameter vector and constraints.params_cand = params_true.copy()params_cand.loc["delta", "value"] = 0.93params_cand["lower_bound"] = [0.79, 0.066, -0.11, 1.04, 0, 0, 0]params_cand["upper_bound"] = [1, 0.072, -0.095, 1.055, 0.5, 0.5, 0.5]constr = constr_base.copy()constr.remove({"loc": "delta", "type": "fixed"})

[65]:

criterion_msm = rp.get_moment_errors_func( params_true, options_new_seed, calc_moments, replace_nans, moments_obs, W)rslt = minimize( criterion=criterion_msm, params=params_cand, algorithm="scipy_lbfgsb", constraints=constr,)rslt

[65]:

{'solution_x': array([0.93]), 'solution_criterion': 7204.902739425344, 'solution_derivative': array([0.]), 'solution_hessian': None, 'n_criterion_evaluations': 1, 'n_derivative_evaluations': None, 'n_iterations': 0, 'success': True, 'reached_convergence_criterion': None, 'message': 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL', 'solution_params': lower_bound upper_bound value category name delta delta 0.790 1.000 0.930 wage_fishing exp_fishing 0.066 0.072 0.070 nonpec_fishing constant -0.110 -0.095 -0.100 nonpec_hammock constant 1.040 1.055 1.046 shocks_sdcorr sd_fishing 0.000 0.500 0.010 sd_hammock 0.000 0.500 0.010 corr_hammock_fishing 0.000 0.500 0.000}

References

• Adda, J., & Cooper, R. W. (2003). Dynamic Economics: Quantitative Methods and Applications. MIT press.

• Adda, J., Dustmann, C., & Stevens, K. (2017). The Career Costs of Children. Journal of Political Economy, 125(2), 293-337.

• Andrews, I., Gentzkow, M., & Shapiro, J. M. (2017). Measuring the Sensitivity of Parameter Estimates to Estimation Moments. The Quarterly Journal of Economics, 132(4), 1553-1592.

• Bruins, M., Duffy, J. A., Keane, M. P., & Smith Jr, A. A. (2018). Generalized Indirect Inference for Discrete Choice Models. Journal of econometrics, 205(1), 177-203.

• Davidson, R., & MacKinnon, J. G. (2004). Econometric Theory and Methods (Vol. 5). New York: Oxford University Press.

• Evans, R. W. (2018, July 5). Simulated Method of Moments (SMM) Estimation. Retrieved November 30, 2019, from https://notes.quantecon.org/submission/5b3db2ceb9eab00015b89f93.

• Frazier, D. T., Oka, T., & Zhu, D. (2019). Indirect Inference with a Non-Smooth Criterion Function. Journal of Econometrics, 212(2), 623-645.

• Gourieroux, M., & Monfort, D. A. (1996). Simulation-based econometric methods. Oxford university press.

• McFadden, D. (1989). A Method of Simulated Moments for Estimation of Discrete Response Models without Numerical Integration. Econometrica: Journal of the Econometric Society, 995-1026.