🔬 Practice: MSM estimation of the retirement model

🚧 🚧 🚧 Not the final version 🚧 🚧 🚧

Today we estimate the retirement model of Class 18 by the method of simulated moments of Class 19: simulated data, moments, the criterion, the disutility of work and the discount factor recovered, and the standard errors checked against a Monte Carlo.

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

What we have and what we write¶

The course architecture again: model object, solver, simulator, estimator. The first three exist, and today’s work is the fourth:

The estimator in msm.py takes the criterion and the simulated moments as function arguments, so the three functions we write are the whole interface between the model and the method.

import sys, time, pathlib, tempfile, urllib.request
import numpy as np
import matplotlib.pyplot as plt
CODE = pathlib.Path('../gh-code/session20-nov3')          # the DC-EGM code of Class 18 and the estimator of Class 19
if not CODE.exists():                                     # otherwise fetch the modules from the code repository
    CODE = pathlib.Path(tempfile.mkdtemp())
    url = 'https://raw.githubusercontent.com/fediskhakov/sb-dse-code/main/session20-nov3/'
    for f in ('retirement.py', 'msm.py'):
        urllib.request.urlretrieve(url + f, CODE / f)
sys.path.insert(0, str(CODE))
from retirement import retirement, solve_dcegm, simulate
from msm import estimate, jacobian, vcov, jtest

---------------------------------------------------------------------------
HTTPError                                 Traceback (most recent call last)
Cell In[1], line 9
      5 if not CODE.exists():                                     # otherwise fetch the modules from the code repository
      6     CODE = pathlib.Path(tempfile.mkdtemp())
      7     url = 'https://raw.githubusercontent.com/fediskhakov/sb-dse-code/main/session20-nov3/'
      8     for f in ('retirement.py', 'msm.py'):
----> 9         urllib.request.urlretrieve(url + f, CODE / f)
     10 sys.path.insert(0, str(CODE))
     11 from retirement import retirement, solve_dcegm, simulate
     12 from msm import estimate, jacobian, vcov, jtest

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

The “data” are simulated as in Class 18, with a seed of their own: 1000 workers with initial wealth uniform on $[2, 20]$, taste shocks $\sigma = 0.1$ and income shocks $\sigma_y = 0.2$, over $T = 10$ periods. The distribution of initial wealth is treated as part of the model, known to the estimator. The parameters we pretend not to know:

def model(theta):
    '''The retirement model at the parameters we estimate; everything else at the values of Class 18'''
    delta, beta = theta
    return retirement(beta=beta, delta=delta, sigma=0.1, sigma_y=0.2, nquad=10)

theta0 = np.array([0.5, 0.95])                            # true delta and beta
N = 1000
rng = np.random.default_rng(3)
M0 = rng.uniform(2, 20, N)                                # initial wealth of the sample
m0 = model(theta0)
Md, Cd, Dd = simulate(m0, solve_dcegm(m0), M0, seed=rng)   # the data: wealth, consumption, work status
print(f'share working by period  {np.round(Dd.mean(0), 2)}')
print(f'mean wealth by period    {np.round(Md.mean(0), 1)}')

Every column is a period of life: the share working falls from 83% to zero in the last period, where no one works, and wealth is run down from 11 to below 2.

Moments¶

The two profiles above are the moments: the share working in periods $0, \dots, T-2$ and mean wealth in periods $1, \dots, T-1$, $q = 18$ numbers. The last period of work status is left out because it is zero whatever the parameters, and the first period of wealth because it is the observed initial condition.

The function returns one row per individual, the contributions $h(z_i)$, so that the data moments are a column mean and their covariance $\Omega$ is one call of np.cov:

def moments(M, C, D):
    '''Moment contributions h(z_i), one row per individual: work status in periods 0..T-2, wealth in periods 1..T-1'''
    return np.column_stack([D[:, :-1], M[:, 1:]])

h = moments(Md, Cd, Dd)                                   # N x q
mhat, Omega = h.mean(0), np.cov(h.T)                      # data moments and their covariance
q, K = mhat.size, theta0.size
Wdiag = np.diag(1 / np.diag(Omega))                       # first-step weighting: inverse variances
print(f'{q} moments, {K} parameters; standard errors of the data moments:')
print(f'  work status  {np.round(np.sqrt(np.diag(Omega)[:9] / N), 3)}')
print(f'  wealth       {np.round(np.sqrt(np.diag(Omega)[9:] / N), 2)}')

The two kinds of moments live on different scales: a residual in a work share is of order 0.01, a residual in mean wealth of order 0.1, and the inverse variances in Wdiag put them on the same footing. With the identity matrix the wealth residuals would decide everything and the work shares nothing.

The simulator with fixed draws¶

The simulated moments solve the model at $\theta$ and simulate $S$ copies of the data, initial wealth included, with the same seed at every $\theta$, as in the criterion of Class 19. The simulated sample is independent of the data, which the asymptotics require.

The parameter walls keep the optimizer where the solver works: $\beta < 1$ and $\delta > 0$.

def simmoments(theta, S=5, seed=1):
    '''Simulated moments at theta: S copies of the data, initial wealth included, same draws at every theta'''
    delta, beta = theta
    if not (0 < delta < 5 and 0.5 < beta < 1): return None   # outside: no solution, or no retirement at all
    m = model(theta)
    sol = solve_dcegm(m)                                  # about 5 milliseconds
    rng = np.random.default_rng(seed)                     # fixed draws: the same generator state at every theta
    Ms, Cs, Ds = simulate(m, sol, rng.uniform(2, 20, S * N), seed=rng)
    return moments(Ms, Cs, Ds).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

t0 = time.perf_counter()
Q0 = criterion(theta0, Wdiag)
print(f'criterion at the truth {Q0:.4f}, one evaluation with S = 5 in {time.perf_counter() - t0:.3f} seconds')

The criterion at the true parameters is not zero: the data moments carry sampling noise and the simulated moments carry simulation noise, and the criterion measures their distance. The figure cuts the criterion along $\delta$ and along $\beta$ through the truth, and zooms in on a window of $\delta$ one thousandth wide.

fig, ax = plt.subplots(1, 3, figsize=(16, 4.5))
dg, bg, dz = np.linspace(0.3, 0.7, 41), np.linspace(0.9, 0.99, 37), np.linspace(0.4995, 0.5005, 51)
for S, col in ((1, 'C0'), (5, 'C1')):
    ax[0].plot(dg, [criterion([d, theta0[1]], Wdiag, S) for d in dg], color=col, lw=1, label=f'S = {S}')
    ax[1].plot(bg, [criterion([theta0[0], b], Wdiag, S) for b in bg], color=col, lw=1)
    ax[2].plot(dz, [criterion([d, theta0[1]], Wdiag, S) for d in dz], '.-', color=col, ms=3, lw=0.7)
ax[0].set(xlabel='δ', yscale='log', title=f'Criterion along δ, β = {theta0[1]}')
ax[1].set(xlabel='β', yscale='log', title=f'Criterion along β, δ = {theta0[0]}')
ax[2].set(xlabel='δ', title='Zoom: the steps of the indicator moments')
ax[0].legend()
for a in ax: a.axvline(theta0[0] if a is not ax[1] else theta0[1], color='k', lw=0.5, ls=':'); a.grid(color='0.9')
plt.show()

• Both cuts have a clear minimum at the truth, three orders of magnitude deep: $\delta$ and $\beta$ are identified by these moments.

• The zoom shows steps: the work status is an indicator $\mathbf 1(\tilde u < P)$, so a simulated share only changes when a fixed draw crosses the moving threshold, by $1/(SN)$ at a time. The steps are smaller with $S = 5$ than with $S = 1$, and they never disappear.

• The steps are why estimate uses Nelder–Mead and why jacobian takes finite differences with a step $h = 10^{-2}$, wide enough to average over many of them: with $h = 10^{-3}$ the derivative of a work share is one step divided by $h$, a number that has nothing to do with the slope.

Estimation and standard errors¶

The first step minimizes the criterion with the diagonal weighting matrix from a deliberately poor start, and then restarts Nelder–Mead from its own solution: on a criterion with steps the simplex can collapse early, and a restart from the solution with a fresh simplex costs a second or two and catches it.

t0 = time.perf_counter()
th1, Q1, nfev = estimate(criterion, x0=(0.3, 0.9), args=(Wdiag,))      # first step, from a poor start
th1, Q1, nfev2 = estimate(criterion, x0=th1, args=(Wdiag,))           # restart from the solution
se1 = np.sqrt(np.diag(vcov(jacobian(simmoments, th1, S=5), Wdiag, Omega, N, S=5)))
print(f'first step   delta {th1[0]:.4f} ({se1[0]:.4f})  beta {th1[1]:.4f} ({se1[1]:.4f})  '
      f'{nfev + nfev2} evaluations in {time.perf_counter() - t0:.1f} s')

The first step lands within a tenth of a standard error of the truth in both parameters. On other data sets drawn from the same model it can land two standard errors away, and that is not a bug either: eighteen moments from an absorbing retirement process and a persistent wealth path are strongly correlated, and a diagonal $W$ treats them as independent, so the estimate follows whichever noise the sample has in common across periods.

The second step uses the optimal weighting matrix $\hat\Omega^{-1}$, starts from the first-step estimate, and delivers the $J$-test with $q - K = 16$ degrees of freedom.

t0 = time.perf_counter()
Wopt = np.linalg.inv(Omega)
th2, Q2, nfev = estimate(criterion, x0=th1, args=(Wopt,))
th2, Q2, nfev2 = estimate(criterion, x0=th2, args=(Wopt,))
se2 = np.sqrt(np.diag(vcov(jacobian(simmoments, th2, S=5), Wopt, Omega, N, S=5)))
J, p = jtest(Q2, N, S=5, df=q - K)
print(f'second step  delta {th2[0]:.4f} ({se2[0]:.4f})  beta {th2[1]:.4f} ({se2[1]:.4f})  '
      f'{nfev + nfev2} evaluations in {time.perf_counter() - t0:.1f} s')
print(f'J-test {J:.1f} on {q - K} degrees of freedom, p-value {p:.3f}')

The second step lands within half a standard error of $\delta = 0.5$ and $\beta = 0.95$, with standard errors three times smaller than the first step for $\delta$ and five times smaller for $\beta$. That is what the optimal $W$ buys when the moments are correlated: it weights the differences between neighboring periods, which the diagonal $W$ cannot see.

The $J$-test does not reject the model that generated the data: $J = 16.9$ on 16 degrees of freedom, $p \approx 0.39$. It does not always come out this way. On nine data sets drawn from the same model with different seeds the p-values ranged from 0.02 to 0.99, with one rejection at the 5% level, and one rejection in twenty is what a 5% test does. Altonji & Segal (1996) is also in play: $\hat\Omega$ has a condition number of $2 \cdot 10^4$ and is estimated from one thousand individuals, and tests based on its inverse over-reject in samples of this size.

The table every MSM paper carries is the fit of the moments, here as a figure: data moments with two standard error bands against the simulated moments at the estimate.

ms = simmoments(th2)
se_m = np.sqrt(np.diag(Omega) / N)
fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
for a, sl, off, lab in ((ax[0], slice(0, 9), 0, 'share working'), (ax[1], slice(9, 18), 1, 'mean wealth')):
    t = np.arange(9) + off
    a.errorbar(t, mhat[sl], yerr=2 * se_m[sl], fmt='o', ms=4, capsize=3, label='data, ± 2 s.e.')
    a.plot(t, ms[sl], 's--', ms=4, lw=0.8, color='C3', label='model at the estimate')
    a.set(xlabel='period', title=f'Fit of the moments: {lab}'); a.grid(color='0.9'); a.legend()
plt.show()

Both profiles are matched within the bands in every period, which is what the $J$-test summarizes. When the test does reject, the figure usually shows why: the optimal $W$ does not close the gaps one by one, it minimizes a quadratic form in which a run of positively correlated gaps on the same side of the data is cheap, and such a run is what the $J$-test picks up.

Monte Carlo check of the standard errors¶

The standard errors come from a formula with a numerical Jacobian computed over steps. The honest check is a Monte Carlo: draw ten data sets from the true model with different seeds, initial wealth included, run the first-step estimator on each, and compare the spread of the estimates to the formula.

Every replication keeps the same simulation seed in simmoments, so only the data change, and the first step is enough for the check because its formula and the second step’s share every ingredient but $W$.

t0 = time.perf_counter()
sol0, ests = solve_dcegm(m0), []
for r in range(10):
    rng_r = np.random.default_rng(100 + r)
    Mr, Cr, Dr = simulate(m0, sol0, rng_r.uniform(2, 20, N), seed=rng_r)   # a new data set from the true model
    hr = moments(Mr, Cr, Dr)
    mh_r, W_r = hr.mean(0), np.diag(1 / np.var(hr, axis=0))
    def crit_r(theta, W):
        '''The criterion with this replication's data moments; the simulation seed stays the same'''
        ms = simmoments(theta)
        return 1e10 if ms is None else (mh_r - ms) @ W @ (mh_r - ms)
    th, _, _ = estimate(crit_r, x0=(0.3, 0.9), args=(W_r,))
    th, _, _ = estimate(crit_r, x0=th, args=(W_r,))
    ests.append(th)
ests = np.array(ests)
print(f'{ests.shape[0]} replications in {time.perf_counter() - t0:.0f} s')
print(f'mean of the estimates      delta {ests[:, 0].mean():.4f}  beta {ests[:, 1].mean():.4f}')
print(f'std of the estimates       delta {ests[:, 0].std():.4f}  beta {ests[:, 1].std():.4f}')
print(f'formula, first step above  delta {se1[0]:.4f}  beta {se1[1]:.4f}')

The estimator is centered on the truth, and the Monte Carlo spread matches the formula within 10% for both parameters, with a Jacobian that averages over the steps. Ten replications is a rough check, not a verdict: the standard error of a standard deviation from ten draws is about a quarter of its value, so agreement within 30% is all the check can promise.

What breaks: the taste shock scale¶

The scale $\sigma$ of the taste shocks was fixed at its true value. Adding it to the parameter vector is one line in model and one in the walls, and it is where the method shows its limits. The criterion along $\sigma$ through the truth, computed with a version of model that takes three parameters:

def model3(theta):
    '''The retirement model with the taste shock scale as a third parameter'''
    delta, beta, sigma = theta
    return retirement(beta=beta, delta=delta, sigma=sigma, sigma_y=0.2, nquad=10)

def criterion3(theta, W, S=5, seed=1):
    if not (0 < theta[0] < 5 and 0.5 < theta[1] < 1 and theta[2] > 0): return 1e10
    m = model3(theta)
    rng = np.random.default_rng(seed)
    Ms, Cs, Ds = simulate(m, solve_dcegm(m), rng.uniform(2, 20, S * N), seed=rng)
    g = mhat - moments(Ms, Cs, Ds).mean(0)
    return g @ W @ g

for s in (0.02, 0.05, 0.08, 0.1, 0.12, 0.15, 0.2):
    print(f'sigma {s:<5} criterion {criterion3([theta0[0], theta0[1], s], Wdiag):.4f}')

Compare with the cuts along $\delta$ and $\beta$ above, which climb more than an order of magnitude within a few hundredths of the truth. Along $\sigma$ the minimum is at the truth, 0.007, but the criterion is flat: it rises to 0.02 at $\sigma = 0.08$ and $\sigma = 0.15$, and to 0.06–0.07 only at $\sigma = 0.05$ and $\sigma = 0.2$. The moments barely respond to $\sigma$, which is the small entry in $D$ that the asymptotics of Class 19 warned about:

• the standard error of $\sigma$ is large, and the finite-difference Jacobian in that direction is dominated by the steps, so even the reported standard error is unreliable

• Nelder–Mead estimating all three parameters stalls on the flat, stepped criterion: started from $\sigma = 0.05$ it stops at $\sigma = 0.068$, from $\sigma = 0.2$ at $\sigma = 0.275$, both with a criterion several times above the minimum, and a restart does not help; only a start near the truth finds $\sigma \approx 0.10$

• on other data sets the sampling noise of eighteen moments can move the flat minimum away from the truth, because the weak signal of $\sigma$ cannot pull it back

What would identify $\sigma$ is a moment that measures how retirement varies given wealth: the dispersion of the retirement age among workers with similar initial wealth, for instance. The practical task asks for it.