Simulation for Inference II — Matching Simulations to Data

36-467/36-667

8 November 2018

\[ \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} \]

In our last episodes

Agenda

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 we’re dealing with well-behaved interior mininma, \[ \hat{\theta}_n \approx \theta^* - (\nabla \nabla M_n(\theta^*))^{-1} \nabla M_n(\theta^*) \] and \[ \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]
logistic.map <- function(x,r) {
  return(4*r*x*(1-x))
}

# 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
  # logistic.map() is defined as above
logistic.map.ts <- 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] = logistic.map(x[t-1],r)
  }
  return(x)
}

Sensitive dependence on initial conditions

traj1 <- logistic.map.ts(1000, r=0.9, initial.cond=0.5001)
traj2 <- logistic.map.ts(1000, 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

par(mfrow=c(1,2))
hist(traj1)
hist(traj2)

logistic.noisy.ts <- function(timelength,r,initial.cond=NULL,noise.sd=0.1) {
    x <- logistic.map.ts(timelength,r,initial.cond)
    return(x+rnorm(timelength,0,noise.sd))
}

mean.logistic <- function(r,n=1e3,s=10) {
    mean(replicate(s,mean(logistic.map.ts(n,r))))
}

mean.logistic.plottable <- function(r,n=1e3,s=10) {
    sapply(r,mean.logistic,n=n,s=s)
}

var.logistic <- function(r,n=1e3,s=10) {
    mean(replicate(s,var(logistic.map.ts(n,r))))
}

var.logistic.plottable <- function(r,n=1e3,s=10) {
    sapply(r,var.logistic,n=n,s=s)
}

Using MSM on the Logistic Map