Housekeeping

• HW12 is now due tomorrow, Friday 4 December, at 6 pm
• HW13 will be due Friday 11 December at 6 pm
• I promise I will in fact grade all the after-class exercises

Model-Based Forecasting

• What?
  • Build a parametric, or at least highly structured, model of the dynamics and the measurement process
    • There’s an internal state \(S(t)\) which evolves according to some model (usually a Markov model)
    • The model for the dynamics of \(S(t)\) is where we try to build in the “physics” (or biology or economics or ecology…)
    • There’re observables \(X(t)\) which are functions of \(S(t)\) plus noise
    • The observation model is where we try to build in our knowledge of the measurement process
    • Both \(S(t)\) and \(X(t)\) might be big vectors, perhaps with spatial structure, etc.
  • Estimate the parameters which show up in the process and observation models
  • Estimate the state based on the history of observations
  • Run the model it forward
  • Run the model forward many times to get at uncertainty
• Why?
  • Non-stationarity and trends
  • Efficient use of information in the data
  • Checking whether a model is right

Fixed Parameters, Fixed State

• Say we know the current state \(S(t)\) and the parameters \(\theta\)
• We can then run the model forward (i.e., simulate it)
  • Get a sequence of future states \(S(t+1), S(t+2), \ldots S(t+T)\)
  • With the states get a sequence of observations \(X(t+1), X(t+2), \ldots X(t+T)\)
• But even with fixed state and parameters, the model will usually produce multiple possible futures
• \(\Rightarrow\) Run the model many times to see the distribution of possible trajectories
  • Just plotting all the \(X\) trajectories can be very informative, but it may also be just a “spaghetti plot”

Fixed Parameters, Unknown State

• The state is often not known exactly
  • What does the monkey have in mind for its hand?
  • How many lynxes and hares are there?
  • How many people are really in compartments \(S\), \(I\), \(R\)?
• We can estimate the state
  • E.g., \(\mathbb{E}\left[ S(t)|X(1), \ldots X(t) \right]\)
• But there’s some uncertainty around this estimate
• Ideally we’d use \(\mathbb{P}\left( S(t)|X(1), \ldots X(t) \right)\)
• Ensemble forecasting: Draw samples from the state distribution and run each one forward
  • This is what weather forecasting does
  • Today’s predicted probability of snow in Pittsburgh on 21 December 2020 = what fraction of runs in the ensemble, started today, have evolved to snow on 21 December 2020?
  • This should remind you of the particle filter…

Learning About the State as We Go

• Start with an ensemble forecast:
  • Many guesses initially about \(S(t)\)
  • Extrapolate forward to guesses about \(S(t+1), S(t+2), \ldots\)
  • Leading to guesses about \(X(t+1), X(t+2), \ldots\), which are the predictions
• Now we actually observe \(X(t+1)\)
• Data assimilation: Use \(X(t+1)\) to update our guesses for \(S(t+1)\) and re-do the forecasting going forward
  • This is exactly what the particle filter does
  • Lots of numerical schemes for data assimilation to try to run it faster at some approximation cost (Koyama et al. 2010)

Unknown Parameters, Unknown State

• The parameters are usually not known exactly
• Different parameters \(\Rightarrow\) different predictions
• How to deal with this?
  1. “Propagation of error”
  2. Bootstrap-type approaches

“Propagation of error”

• If we have a parameter estimate \(\hat{\theta}\) and standard error \(\hat{se}\), then estimate \(f(\theta)\) as \(f(\hat{\theta})\)
• An approximate standard error for \(f(\theta)\) would be \(\left|\frac{d f}{d\theta}\right| \hat{se}\)
  • Approximation becomes exact as \(\hat{se} \rightarrow 0\)
• Extends to multiple arguments: we estimate \(f(\theta_1, \theta_2, \ldots \theta_p)\) by \(f(\hat{\theta}_1, \hat{\theta}_2, \ldots \hat{\theta}_p)\), with standard error \[ \approx \sqrt{\sum_{i=1}^{p}{\left(\frac{\partial f}{\partial \theta_i}\right)^2 \hat{se}_i^2} + 2\sum_{i=1}^{p}{\sum_{j=i+1}^{p-1}{\frac{\partial f}{\partial \theta_i}\frac{\partial f}{\partial \theta_j} \mathrm{Cov}[\hat{\theta}_i, \hat{\theta}_j]}}} \]
  • Just a Taylor expansion; see handout on website
• Lets you propagate uncertainty from parameters to functions of parameters (like predictions) just by taking derivatives
  • Taking derivatives in a big, complicated model can be hard!
• If you’re willing to believe \(\hat{\theta}\) has a Gaussian distribution, then \(f(\hat{\theta})\) does too (approximately)
  • This is called the delta method

Bootstrap approaches to unknown parameters

• Bootstrap the training data
  • Model-based bootstraps recommended when dealing with non-stationary data
• Estimate parameters on each bootstrap run
• For each parameter estimate, run your forecast as if those parameter values were fixed
  • May want to run an ensemble forecast for each parameter value
• Aggregate over all simulation runs

Model Checking and Diagnostics

• Before relying on the model’s predictions, check the model
  • Usually a good idea before wasting time on things like confidence sets of parameters
• You remember a lot of “diagnostics” for linear regression models from your regression class
• For dynamic models, the basic diagnostic is that \[ \frac{X(t) - \mathbb{E}\left[ X(t)|X(1), \ldots X(t-1) \right]}{\sqrt{\mathrm{Var}\left[ X(t)|X(1), \ldots X(t-1) \right]}} \] should be mean 0 with variance 1
• If \(X(t)\) always has a continuous distribution given the past, say with CDF \(F_t(r) \equiv \mathbb{P}\left( X(t) \leq r| X(1), \ldots X(t-1) \right)\), then the sequence \(Q(t) = F_t(X(t))\) should be uniformly distributed on \([0,1]\) and IID
  • Similar but more complicated rules for discrete or discontinuous distributions
• The basic idea: once you’ve removed or adjusted for everything the model predicts, what’s left should be pure noise

Model Checking by Comparisons

• There may be multiple models (perhaps variations on a theme)
• One way to chose would be a pure prediction contest, i.e., try to get the best forecast of a common set of variables
• May be more illuminating to look for qualitative distinctions which you can compare to the data (are there cycles? a hard lower or upper bound on some process?)
• May be even more illuminating to look for new phenomena/variables
• This gets us in to philosophy of science (Glymour 1980; Mayo 1996) but sometimes that can’t be helped
• Alternatively: present results from all the models which you can’t rule out

Summing Up

1. With fixed state and fixed parameters, run the model forward multiple times to get multiple simulated trajectories
2. With unknown state but fixed parameters, get an ensemble of states and go to (1)
   • Updating the state as we see more = data assimilation, often through the particle filter
3. With unknown state and unknown parameters, get an ensemble of parameters and go to (2)
4. Try to do model checking before investing lots of time into model-based forecasting

Backup: Visualizing Many Simulation Runs

• Every run of the simulation gives us a trajectory, which is a function of time; that’s why I write it \(X(t)\)
• The easiest way to display all these at once is just to plot them one atop the other
• We’ve been doing this without remarking on it, as with the SIR model we looked at in Lecture 20

Backup: Visualizing Many Simulation Runs

• With many runs, this can get kind of messy (“spaghetti plots”):

• There are many options for summarizing things like this in more comprehensible and visually-pleasing ways. Common ones:
  1. Mean \(\pm 2\) standard deviations, computed at each value of \(t\)
  2. Median plus high and low quantiles, computed at each value of \(t\)
  3. “Functional boxplots”

Backup: Visualizing Many Simulation Runs

Mean \(\pm 2\) standard deviations, at each value of \(t\)

Comments:

• Fine idea if everything’s Gaussian
• Still informative in many other circumstances but may also give un-“physical” (or un-“biological”) limits, like the negative values here
  • OTOH this does tell us the variance at those times is really big

Backup: Visualizing Many Simulation Runs

Median and high and low quantiles, at each value of \(t\)

• Comments:
  • Common to pick say 25% and 75% quantiles, getting a 50% probability interval
  • Or 5% and 95% to get a 90% probability interval
  • Or, as in the plot, both (using shading or color to reinforce the sense of moving away from the median)
• Drawbacks of both the time-by-time quantile and time-by-time mean approach:
  • The typical trajectory isn’t, necessarily, near the typical value at all times!
  • These “fixed-time” summaries can give a misleading impression of the size and timing of extremes, if the timing of extremes is itself variable
  • E.g., in an epidemic model like this, where the time of peak infection varies, the fixed-time quantiles can miss the peak (Juul et al. 2020)

Backup: Visualizing Many Simulation Runs

“Functional boxplots” (Sun and Genton 2011) or “curve boxplots” (Mirzargar, Whitaker, and Kirby 2014)

• Look at the distance between entire curves
• Try to find the “deepest” curve, the one most nearly at the center of the distribution
• Order curves outward by depth from that deepest curve
• Stop when you have, say, 90% of the curves
• Draw the resulting tube around the deepest curve, and highlight the outliers falling beyond it (as in a regular boxplot)