📖 Method of simulated moments

🚧 🚧 🚧 Not the final version 🚧 🚧 🚧

EGM solves the consumption-savings model in milliseconds, but there is no likelihood for consumption and wealth data. The method of simulated moments (MSM) needs only the simulator: pick a few statistics of the data, and find the parameters at which the model reproduces them. Theory today, practice on Tuesday.

cd sb-dse-code
git reset --hard HEAD
git pull

From solver and simulator to estimator¶

A solver and a simulator that are debugged let us study a theory: run the model for a range of parameter values and see what the decision maker does. The ultimate goal is to match the model to the observed data, in order to

• quantify the effects of the parts of the theoretical setup

• run counterfactual experiments: simulate behavior under a policy regime that has never been observed

• support or falsify the theory by its fit to the data

The workflow of a structural project is the arc of this course: model, implementation, solver, understanding how the model works, estimation, validation, counterfactuals. Step four is where the time goes: solve, simulate, ask whether the result makes economic sense, repeat. It is also where the moments for step five are found.

Estimation versus calibration¶

The terms are often used interchangeably, and calibration exercises also search for the parameters that fit the data best. What estimation adds:

• standard errors, a measure of the variability of the estimates

• a study of identification, to make sure that no other parameter values explain the data equally well

Estimation often concerns only a subset of the parameters, with the rest fixed at values from the literature, like a calibration. The NFXP estimation of Class 10 fixed $\beta$, for instance.

Why not maximum likelihood?¶

The MLE of Class 10 is the most efficient estimator there is, when the model gives the probability of the observed data. In the bus engine model it does: the logit choice probabilities. In the consumption-savings model it does not, or not easily:

• consumption is a deterministic function of the state, $c = c(M)$. Observed consumption that is off the policy function has likelihood zero, so a likelihood needs measurement error in consumption and an assumption about its distribution

• wealth is a sum of past shocks, so the likelihood of a wealth path is an integral over every unobserved income draw, and in a panel of $T$ periods that is a $T$-dimensional integral per individual

• the data are often aggregated: wealth by age from one survey, consumption from another, and no individual with both

Simulation is cheap where integration is not. The simulator of Class 16 produces a synthetic panel for any $\theta$, and any statistic of the synthetic panel can be compared to the same statistic of the data.

Method of moments, GMM and simulated moments¶

Matching statistics of the model to statistics of the data has a long history:

• method of moments: as many moments as parameters, a system of equations, solved exactly

• generalized method of moments (GMM) of Hansen (1982): more moments than parameters, minimize a weighted distance between the data moments and the theoretical moments

• method of simulated moments (MSM) of McFadden (1989): GMM with the theoretical moments replaced by averages over simulated data

Let $z_i$, $i = 1, \dots, N$, be the observed data, one record per individual, and let $h(z)$ be a vector of $q$ statistics of one record, the moment functions. Three objects share the name “moments”:

The model moments are the expectation of $h$ under the model with parameters $\theta$, the object GMM works with. MSM replaces them by the average over a synthetic sample $\tilde z_s(\theta)$ of $S N$ records, simulated from the model at $\theta$ with random draws $\tilde u_s$. The vector

of $q$ moment conditions is zero in expectation at the true $\theta_0$ and at no other point, if the moments identify the parameters.

The last clause is the whole trick. With the draws fixed, every evaluation of $Q_N(\theta)$ solves and simulates the model on the same shocks, so the simulated moments are a deterministic function of $\theta$ that an optimizer can minimize.

Redraw the shocks at every $\theta$ and the criterion is a random function: an optimizer cannot tell a change in the parameters from a change in the draws.

Three conditions make the estimator work, in the words of Corbae & Guvenen (2024):

• the simulation draws are independent of the data

• the model is identified: $m(\theta) = m(\theta_0)$ only at $\theta_0$

• the simulated moments are continuous in the mean: $\tilde m_S(\theta) \to m(\theta)$ as the simulated sample grows

The count of moments against parameters is the order condition, necessary but not sufficient:

• $q > K$: overidentified, the usual case, with a test of the model for free

• $q = K$: just identified, the criterion is zero at the optimum and the weighting matrix does not matter

• $q < K$: underidentified, no estimator

The simulated criterion¶

The data are simulated too, from the “true” parameters, so that we know the answer. The module deaton_egm.py holds the model object of Class 16, the EGM solver of Class 17, which returns the policy as points on the endogenous grid, and the simulator.

import sys, time, pathlib, tempfile, urllib.request
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.stats import chi2
CODE = pathlib.Path('../gh-code/session19-oct29')         # the model and EGM solver of Classes 16-17
if not CODE.exists():                                     # otherwise fetch it from the code repository
    CODE = pathlib.Path(tempfile.mkdtemp())
    urllib.request.urlretrieve('https://raw.githubusercontent.com/fediskhakov/sb-dse-code/main/session19-oct29/deaton_egm.py',
                               CODE / 'deaton_egm.py')
sys.path.insert(0, str(CODE))
from deaton_egm import deaton, solve_egm, simulate

---------------------------------------------------------------------------
HTTPError                                 Traceback (most recent call last)
Cell In[1], line 9
      5 from scipy.stats import chi2
      6 CODE = pathlib.Path('../gh-code/session19-oct29')         # the model and EGM solver of Classes 16-17
      7 if not CODE.exists():                                     # otherwise fetch it from the code repository
      8     CODE = pathlib.Path(tempfile.mkdtemp())
----> 9     urllib.request.urlretrieve('https://raw.githubusercontent.com/fediskhakov/sb-dse-code/main/session19-oct29/deaton_egm.py',
     10                                CODE / 'deaton_egm.py')
     11 sys.path.insert(0, str(CODE))
     12 from deaton_egm import deaton, solve_egm, simulate

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:240, in urlretrieve(url, filename, reporthook, data)
    223 """
    224 Retrieve a URL into a temporary location on disk.
    225 
   (...)    236 data file as well as the resulting HTTPMessage object.
    237 """
    238 url_type, path = _splittype(url)
--> 240 with contextlib.closing(urlopen(url, data)) as fp:
    241     headers = fp.info()
    243     # Just return the local path and the "headers" for file://
    244     # URLs. No sense in performing a copy unless requested.

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:215, in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    213 else:
    214     opener = _opener
--> 215 return opener.open(url, data, timeout)

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:521, in OpenerDirector.open(self, fullurl, data, timeout)
    519 for processor in self.process_response.get(protocol, []):
    520     meth = getattr(processor, meth_name)
--> 521     response = meth(req, response)
    523 return response

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:630, in HTTPErrorProcessor.http_response(self, request, response)
    627 # According to RFC 2616, "2xx" code indicates that the client's
    628 # request was successfully received, understood, and accepted.
    629 if not (200 <= code < 300):
--> 630     response = self.parent.error(
    631         'http', request, response, code, msg, hdrs)
    633 return response

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:559, in OpenerDirector.error(self, proto, *args)
    557 if http_err:
    558     args = (dict, 'default', 'http_error_default') + orig_args
--> 559     return self._call_chain(*args)

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:492, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
    490 for handler in handlers:
    491     func = getattr(handler, meth_name)
--> 492     result = func(*args)
    493     if result is not None:
    494         return result

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/urllib/request.py:639, in HTTPDefaultErrorHandler.http_error_default(self, req, fp, code, msg, hdrs)
    638 def http_error_default(self, req, fp, code, msg, hdrs):
--> 639     raise HTTPError(req.full_url, code, msg, hdrs, fp)

HTTPError: HTTP Error 404: Not Found

theta0 = np.array([0.92, 0.5])                            # true beta and sigma, to be recovered
N, T, tail = 500, 40, 10                                  # individuals, simulated periods, observed periods
m0 = deaton(beta=theta0[0], sigma=theta0[1])              # R = 1.05, 100 grid points on [0, 20]
rng = np.random.default_rng(2026)
M0 = np.exp(rng.standard_normal(N))                       # initial wealth, log-normal
Wd, Cd = simulate(m0, solve_egm(m0), M0, T, seed=rng)     # the "data": N wealth and consumption paths

def moments(W, C):
    '''Moment functions h(z), one row per individual, from the last `tail` periods of the paths'''
    Wt, Ct = W[:, -tail:], C[:, -tail:]                   # the first 30 periods are burn-in, not observed
    return np.column_stack([Wt.mean(1), Wt.std(1), np.diff(np.log(Ct), axis=1).std(1)])

h = moments(Wd, Cd)                                       # N x q individual contributions
mhat = h.mean(0)                                          # data moments
Omega = np.cov(h.T)                                       # covariance of the contributions, q x q
for name, mh, s in zip(('mean wealth', 'std wealth', 'std consumption growth'), mhat, np.sqrt(np.diag(Omega) / N)):
    print(f'{name:<25} {mh:7.4f}  (s.e. {s:.4f})')

Only the last ten of the forty simulated periods are “observed”: the first thirty wash out the initial wealth, so the data come from the ergodic distribution of the model. The simulation below runs the same forty periods.

The covariance matrix $\Omega$ of the individual contributions $h(z_i)$ is the one object the asymptotics need from the data, and it comes for free once the moments are means over individuals.

The simulated moments solve the model at $\theta$ and simulate $S$ copies of the data with a fixed seed, so the same draws are used at every $\theta$:

def simmoments(theta, S=5, seed=1):
    '''Simulated moments at theta: S copies of the data, same draws at every theta'''
    beta, sigma = theta
    if not (0 < beta < 1 / 1.05 and sigma > 0): return None   # EGM does not converge for beta R >= 1
    m = deaton(beta=beta, sigma=sigma)
    pol = solve_egm(m)                                    # about 3 milliseconds
    Ws, Cs = simulate(m, pol, np.tile(M0, S), T, seed=np.random.default_rng(seed))   # fixed draws
    return moments(Ws, Cs).mean(0)

def criterion(theta, W, S=5, seed=1):
    '''The MSM criterion g' W g'''
    ms = simmoments(theta, S, seed)
    if ms is None: return 1e10                            # a wall for the optimizer
    g = mhat - ms
    return g @ W @ g

The figure plots the criterion in $\beta$, with $\sigma$ at its true value and a diagonal weighting matrix, for three simulators: $S = 1$ with fixed draws, $S = 10$ with fixed draws, and $S = 1$ with fresh draws at every $\beta$.

Wdiag = np.diag(1 / np.diag(Omega))                       # moments scaled by their variances
fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
for a, (lo, hi) in zip(ax, ((0.85, 0.95), (0.912, 0.928))):
    bg = np.linspace(lo, hi, 81)
    a.plot(bg, [criterion([b, theta0[1]], Wdiag, S=1) for b in bg], lw=1, label='S = 1, fixed draws')
    a.plot(bg, [criterion([b, theta0[1]], Wdiag, S=10) for b in bg], lw=1, label='S = 10, fixed draws')
    a.plot(bg, [criterion([b, theta0[1]], Wdiag, S=1, seed=k) for k, b in enumerate(bg)], '.', ms=3, label='S = 1, fresh draws')
    a.axvline(theta0[0], color='k', lw=0.5, ls=':')
    a.set(xlabel='β', title=f'MSM criterion in β on [{lo}, {hi}], σ = {theta0[1]}'); a.grid(color='0.9')
ax[0].set(yscale='log'); ax[0].legend()
plt.show()

Three things to see:

• with fixed draws the criterion is a smooth function of $\beta$ and has one minimum, close to the true value

• with fresh draws it is a cloud: in the zoom the simulation noise is as large as the differences the optimizer needs to see, and the minimum is wherever the draws happened to fall

• a larger $S$ lowers the noise of the simulated moments, so the fixed-draw criterion with $S = 10$ is centered better than with $S = 1$ — but it is the fixing, not the number of draws, that makes the function minimizable

Where it breaks: a discrete choice in the simulator. The retirement decision of Class 18 is an indicator $\mathbf 1(\tilde u < P)$, so a simulated retirement rate is a step function of $\theta$ even with fixed draws: it only changes when a draw crosses the threshold.

This is the problem of McFadden (1989), and the reason the asymptotics of Pakes & Pollard (1989) do not assume a smooth criterion. Gradient-based optimizers fail on the steps; derivative-free ones, Nelder–Mead in our case, walk over them as long as the steps are small compared to the simplex.

Asymptotic properties¶

The MSM estimator is a GMM estimator with a noisy version of the model moments, and it inherits the GMM theory of Hansen (1982) with one correction. Let $D$ be the $q \times K$ Jacobian of the model moments,

and $\Omega$ the $q \times q$ asymptotic covariance of the data moments, $\sqrt N (\hat m_N - m(\theta_0)) \xrightarrow{d} N(0, \Omega)$.

Where does the factor $1 + 1/S$ come from? The moment conditions at $\theta_0$ carry two independent sources of noise,

the data moments with variance $\Omega / N$ and the simulated moments with variance $\Omega / (SN)$, because they average $SN$ records from the same distribution. Their sum has variance $(1 + 1/S)\, \Omega / N$, and this is the only place where the simulation enters the theory.

What the factor costs: with $S = 1$ the standard errors are $\sqrt 2$ times the GMM ones, with $S = 10$ they are 5% larger, and simulating more than ten copies of the data buys very little.

Three more things to read off the theorem:

• Consistency needs no good $W$ and no large $S$. The estimator is consistent with the identity matrix and $S = 1$. What $W$ and $S$ buy is precision.

• The sandwich $(D^T W D)^{-1} D^T W \Omega W D (D^T W D)^{-1}$ is the robust form of Class 10, with the Jacobian of the moments in place of the Hessian of the likelihood. It collapses to $(D^T \Omega^{-1} D)^{-1}$ only at the optimal weighting matrix.

• Local identification is a rank condition: $D$ has to have $K$ linearly independent columns, so every parameter must move at least one moment, and no two parameters may move the moments in the same direction. A flat direction of the criterion is a column of $D$ that is zero or a combination of the others.

• Precision is governed by $D$. A moment that barely responds to a parameter is a small entry of $D$ and a large entry of $(D^T W D)^{-1}$: this is how local identification is linked to the standard errors, and a weakly identified parameter announces itself through them.

Where it breaks: every formula above needs $D$, the derivative of a function we only know through a simulator. The theorem does not require a smooth criterion, Pakes & Pollard (1989) cover step functions, but the standard errors do need the slope.

In practice $D$ is computed by finite differences of $\tilde m_S(\theta)$ with fixed draws, and the step size must be small enough to be a derivative but large enough to jump over the steps and the interpolation kinks of the solver. Tuesday’s class shows how much rides on that choice.

The weighting matrix and the two-step estimator¶

The optimal weighting matrix is the inverse of the covariance of the moment conditions, which is unknown but estimable. For moments that are means over individuals, $\hat\Omega$ is the sample covariance of the contributions $h(z_i)$, the one line np.cov(h.T) above, computed from the data once.

For moments of other kinds, autocovariances of a time series for instance, a Newey–West estimator or a bootstrap replaces it, and when the data moments come from published tables rather than micro data, $\hat\Omega$ is itself simulated from the model at a first-step estimate.

Iterating steps 2 and 3 until $\theta$ stops changing is the iterated MSM; the continuously updated version makes $W$ a function of $\theta$ inside the criterion. In practice the weighting matrix is often simplified, and Altonji & Segal (1996) is the reason: in small samples the optimal $W$ is estimated with error that is correlated with the moment conditions, and the second step can be more biased than the first. The choices seen in applied work, from the simplest up:

• the identity, in the first step of a multi-step estimation

• a diagonal matrix of inverse variances, which brings the moments to a common scale and ignores their covariances

• manually chosen weights, to bring all moments to the same scale or to emphasize the moments the study is about

• the optimal $\hat\Omega^{-1}$, from the data or from a first-step simulation

• an iteratively updated $W$, a Newey–West estimate for serially correlated moments, and model-specific adjustments

Many ways to skin a cat, and the choice should be reported with the estimates.

The test of overidentifying restrictions¶

With more moments than parameters the criterion is not zero at the optimum, and its size says whether the extra moments are compatible with the model.

The test is valid only with the optimal weighting matrix: with any other $W$ the criterion at the optimum has a different distribution. And it tests the moments chosen, not the model: a model that fits three moments can be rejected by a fourth.

def estimate(W, S=5, x0=(0.9, 0.4), seed=1):
    '''Minimize the MSM criterion by Nelder-Mead: derivative free, so steps in the simulator do not matter'''
    res = minimize(criterion, x0, args=(W, S, seed), method='Nelder-Mead',
                   options={'xatol': 1e-5, 'fatol': 1e-9})
    assert res.success, res.message
    return res.x, res.fun, res.nfev

def jacobian(theta, S=5, h=1e-4):
    '''Jacobian D of the model moments by central differences of the simulated moments, q x K'''
    theta, D = np.asarray(theta, float), []
    for e in np.eye(theta.size) * h:
        D.append((simmoments(theta + e, S) - simmoments(theta - e, S)) / (2 * h))   # same draws on both sides
    return np.array(D).T

def vcov(theta, W, S=5):
    '''Asymptotic variance of the MSM estimator, the sandwich with the (1 + 1/S) correction'''
    D = jacobian(theta, S)
    assert np.linalg.matrix_rank(D) == theta.size          # local identification: full column rank
    bread = np.linalg.inv(D.T @ W @ D)
    return (1 + 1 / S) * bread @ D.T @ W @ Omega @ W @ D @ bread / N

S = 5
for name, W in (('identity', np.eye(3)), ('diagonal', Wdiag), ('optimal', np.linalg.inv(Omega))):
    t0 = time.perf_counter()
    th, Q, nfev = estimate(W, S)
    se = np.sqrt(np.diag(vcov(th, W, S)))
    line = f'{name:<9} beta {th[0]:.4f} ({se[0]:.4f})  sigma {th[1]:.4f} ({se[1]:.4f})  {nfev} evaluations, {time.perf_counter() - t0:.2f} s'
    if name == 'optimal':                                 # the J-test needs the optimal weighting matrix
        J = N / (1 + 1 / S) * Q
        line += f'  J = {J:.2f}, p-value {chi2.sf(J, 3 - 2):.3f}'
    print(line)

All three weighting matrices recover $\beta = 0.92$ and $\sigma = 0.5$ within one standard error, in under a second each, and the standard errors shrink as the weighting improves. The $J$-test does not reject: the third moment is compatible with a model that fits the first two, as it should be when the data come from the model itself.

The formula is not the only way to a standard error. A bootstrap resamples the individuals, re-runs the whole two-step estimation with a fresh simulation seed, and reports the spread of the estimates across replications.

It costs hundreds of estimations, minutes here and days in a large model, but it needs no Jacobian and no asymptotic argument. The Monte Carlo of the next class is the same idea with the true parameters known.

Choosing the moments¶

The choice of moments is the crucial part of an MSM estimation, and more art than science. The principles:

• the moments must be sensitive to the parameters: understanding how the model works means understanding what variation the parameters induce in the simulated data, and the $q$ moments should represent that variation

• every parameter needs its own moment, in the sense of the rank condition: the columns of $D$ are the responses of the moments to the parameters, and they must be linearly independent

• moments must be measurable precisely in the data: a moment with a large variance in $\hat\Omega$ gets a small weight and contributes little

• moments should be interpretable, so that the fit can be shown in a table of data moments against simulated moments, which every MSM paper carries

The figure maps the criterion over $(\beta, \sigma)$ with the diagonal weighting matrix, with all three moments on the left and with the two wealth moments only on the right.

bg, sg = np.linspace(0.88, 0.95, 36), np.linspace(0.3, 0.7, 33)
B, Sg = np.meshgrid(bg, sg)
fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
for a, keep, title in zip(ax, (slice(0, 3), slice(0, 2)), ('three moments', 'two wealth moments only')):
    Wk = np.zeros((3, 3)); Wk[keep, keep] = Wdiag[keep, keep]   # drop a moment by zeroing its weight
    Q = np.array([[criterion([b, s], Wk, S=5) for b in bg] for s in sg])
    cs = a.contourf(B, Sg, np.log10(Q), levels=20, cmap='viridis')
    a.plot(*theta0, 'r+', ms=12, mew=2)
    a.set(xlabel='β', ylabel='σ', title=f'log10 of the MSM criterion, {title}')
    fig.colorbar(cs, ax=a)
plt.show()

With three moments the criterion has one basin around the truth. With the two wealth moments only, the basin stretches into a valley: a lower $\beta$ with a higher $\sigma$ produces nearly the same mean and spread of wealth, impatience traded against precautionary saving.

The rank condition still holds along the valley, $D$ has rank two, but its columns are nearly collinear: the standard errors explode and an optimizer stops anywhere along the floor. The consumption moment breaks the trade-off because $\sigma$ moves consumption growth and $\beta$ does not.

Where it breaks, in applications:

• a parameter that moves no chosen moment: a flat criterion in that direction, an optimizer that stops at the starting value, and standard errors from a singular matrix

• two parameters with the same signature on the moments: a valley, as on the right

• moments that are step functions of $\theta$ in the simulator, which can hide a valley from a derivative-free optimizer and ruin the finite-difference Jacobian

Andrews et al. (2017) give a measure of what each moment contributes to each estimate, $\Lambda = -(D^T W D)^{-1} D^T W$, the matrix that maps a perturbation of the moments into a perturbation of the estimates. It is a by-product of the standard errors and belongs in every MSM paper.

MSM in practice¶

Why is the method so widely used in applied work?

• it is not data hungry: it matches aggregated moments, and can combine different sources of data, a wealth profile from one survey and a retirement hazard from another

• it does not rely on distributional assumptions as much as MLE, which needs the whole density of the data

• it only needs the simulator, which the structural project has anyway

The price is efficiency: the standard errors are larger than those of the MLE, the weighting matrix is often simplified because of its small-sample bias, and the choice of moments is a judgment call that has to be defended. Two relatives extend the idea: indirect inference Gourieroux et al., 1993Smith, 1993 matches the coefficients of an auxiliary model, a regression or a VAR, estimated on the data and on the simulations, and simulated moments of time series Lee & Ingram, 1991Duffie & Singleton, 1993 match autocovariances of one long simulated path.

Two life-cycle applications of the kind we estimate next class:

• Gourinchas & Parker (2002) estimate the discount factor and risk aversion of a consumption-savings model from the age profile of consumption, with the age profile of wealth as a check

• French (2005) estimates a labor supply and retirement model from profiles of hours, participation and assets by age, the model of Class 20 with health and wages added

Iskhakov & Keane (2021) estimate a dynamic model of labor supply, savings and human capital on Australian survey data by MSM, and use it to evaluate a reform of the pension system: the whole arc from model to counterfactual in one paper.