Simulation for Inference II — Matching Simulations to Data


27 October 2020 (Lecture 17)

\[ \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Prob}[1]{\mathbb{P}\left[ #1 \right]} \newcommand{\Probwrt}[2]{\mathbb{P}_{#1}\left( #2 \right)} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}`left[ #1 \right]} \newcommand{\Expectwrt}[2]{\mathbb{E}_{#1}\left[ #2 \right]} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator*{\argmax}{argmax} \]

May their memory be a blessing

In our last episodes


How we estimate, in general

\[ \hat{\theta}_n = \argmin_{\theta}{M_n(\theta)} \]

If \(M_n(\theta) \rightarrow m(\theta)\) as \(n\rightarrow\infty\)

and \(m(\theta)\) has a unique minimum at the true \(\theta^*\)

then (generally) \(\hat{\theta}_n \rightarrow \theta^*\) (consistency)

If \(m\) has a well-behaved interior minimum, \[ \hat{\theta}_n \approx \theta^* - (\nabla \nabla M_n(\theta^*))^{-1} \nabla M_n(\theta^*) \] and so \[ \Var{\hat{\theta}_n} \approx (\nabla\nabla m(\theta^*))^{-1} \Var{\nabla M_n(\theta^*)} (\nabla\nabla m(\theta^*))^{-1} \]

The Progress of Statistical Methods

  1. Calculate likelihood, solve explicitly for MLE
  2. Can’t solve for MLE but can still write down likelihood, calculate it, and maximize numerically
  3. Even calculating the likelihood is intractable

Why Finding the Likelihood Becomes Hard

Method of simulated moments

Assume: the data really came from our simulation model, at some true value \(\theta^*\) of the adjustable parameters \(\theta\)

  1. Pick your favorite \(q\)-dimensional vector of statistics \(B\) (“generalized moments”)
  2. Calculate from data, \(B(obs) = B(X)\)
  3. Pick a staring parameter value \(\theta\)
    1. simulate multiple times, say \(s\), getting \(\tilde{X}_1, \ldots \tilde{X}_s\)
    2. calculate average of \(B\), \(\overline{B}(\theta, n, s) \equiv \frac{1}{s}\sum_{i=1}^{s}{B(\tilde{X}_i)}\)
    3. For large \(n\), \(\overline{B}(\theta, n, s) \approx \Expectwrt{\theta}{B}\)
  4. Adjust \(\theta\) so expectations are close to \(B(obs)\)

Method of simulated moments

Method of simulated moments

A challenging example

The logistic map in R

# Do one step of the logistic map, $x(t+1) = 4rx(t)(1-x(t))$
# Inputs: initial condition (x)
  # logistic parameter (r)
# Output: next value of the logistic map
# Presumes: x is a single numerical value in [0,1]
  # r is a single numerical value in [0,1] <- function(x,r) {

# Generate a time series from the logistic map, $x(t+1)=4rx(t)(1-x(t))$
# Inputs: length of desired time series (timelength)
  # logistic parameter (r)
  # value of x(1) (initial.cond, uniformly randomly generated if omitted)
# Output: vector of length timelength
# Presumes: r is a single numerical value in [0,1]
  # timelength is a positive integer
  # initialcond is a single numerical value in [0,1], or omitted
  # is defined as above <- function(timelength, r, initial.cond=NULL) {
  x <-vector(mode="numeric", length=timelength)
  if(is.null(initial.cond)) {
    x[1] <- runif(1)
  } else {
    x[1] <- initial.cond
  for (t in 2:timelength) {
    x[t] =[t-1],r)

Sensitive dependence on initial conditions

traj1 <-, r=0.9, initial.cond=0.5001)
traj2 <-, r=0.9, initial.cond=0.4998)
plot(1:100, traj1[1:100], ylim=c(0,1), xlab="t", ylab="X(t)", type="b")
points(1:100, traj2[1:100], pch=2, col="blue", type="b")

Sensitive dependence on initial conditions

plot(traj1, traj2)
rug(x=traj1, side=1, col="black", ticksize=-0.01) # negative sizes for pointing outward
rug(x=traj2, side=2, col="blue", ticksize=-0.01)

Long-run stability