# Impact of the model error in numerical model prediction

I'm having trouble wrapping my head around characterizing model prediction uncertainties in the case of a calibrated numerical model. Let's assume we are trying to calibrate a set of input parameters $X$ for our model $f$ based on some data $Y$, i.e.:

$Y=f(X)+ε_y$

If we assume a Bayesian viewpoint, built into this equation is one set of fixed variables (our data: $Y$), and two random variables (uncertainty in our parameters: $X$ and the model error: $ε_y$). An MCMC analysis could then be conducted to obtain the uncertainty in $X$ and $ε_y$.

The issue I’m trying to figure out is how do we then use this to make predictions from our numerical model. For example, say we want to generate some future random realization for $Y$. To do this, we need to take into consideration both the uncertainty in the input parameters ($X$) as well as the estimated model errors ($ε_y$). This is relatively straightforward when estimating unseen $Y$ outputs, but becomes more difficult if we want to modify our models to estimate a future condition $Z$.

In numerical modelling terms, what I am trying to get at here is that we calibrate our models ($f$) to some data ($Y$), then after calibration we modify our models ($f→g$) to include the future conditions and want to predict the impact of such conditions ($Z$). The confusion I'm having is that I don’t think we can just take our current model ($f$), modify it by adding the future conditions ($f→g$), and then drive the estimated uncertainty in our parameter set ($X$) through it, i.e.:

$Z=g(X)$

The issue I see with this is that it will cause our uncertainty estimates to be over-fit, as if fails to drive through the full uncertainty in the problem by dropping the modelling error ($ε_y$) term. In reality should the equation not have an additional model error ($ε_z$) with respect to $Z$, i.e.:

$Z=g(X)+ε_z$

However, since it is impossible to infer the error of future unobservable states ($ε_z$), isn’t this problem inherently unknowable?

• Model errors are typically systematic (not random). Your first have to be clear on the distributional properties of $\epsilon_y$ before making predictions with it.... – Pascal Apr 9 '18 at 11:17

## 1 Answer

I'm not sure I'm addressing all of your concerns exactly, but perhaps this will provide a framework to get you moving in the right direction.

Suppose \begin{equation} y_t = f(x_t,\theta) + \varepsilon_t \end{equation} where \begin{equation} \varepsilon_t \stackrel{\text{iid}}{\sim} \textsf{N}(0,\sigma^2) . \end{equation} In this setup, $x_t$ is an observed covariate (a "condition"?) and $\theta$ is an unobserved parameter vector. The observation $y_t$ is normally distributed with mean $f(x_t,\theta)$ and variance $\sigma^2$.

Given $T$ observations (and associated covariates/conditions), we have \begin{equation} p(y_{1:T}|x_{1:T},\theta,\sigma^2) = \prod_{t=1}^T p(y_t|x_t,\theta,\sigma^2), \end{equation} where \begin{equation} p(y_t|x_t,\theta,\sigma^2) = \textsf{N}(y_t|f(x_t,\theta),\sigma^2) . \end{equation} Given a prior distribution for the unobserved parameter, $p(\theta,\sigma^2)$, we can compute the posterior distribution \begin{equation} p(\theta,\sigma^2|y_{1:T},x_{1:T}) = \frac{p(y_{1:T}|x_{1:T},\theta,\sigma^2)\,p(\theta,\sigma^2)}{p(y_{1:T}|x_{1:T})} , \end{equation} where \begin{equation} p(y_{1:T}|x_{1:T}) = \iint p(y_{1:T}|x_{1:T},\theta,\sigma^2)\,p(\theta,\sigma^2)\,d\theta\,d\sigma^2 . \end{equation}

We now turn to prediction. We wish to predict $y_{T+1}$. We will need $x_{T+1}$. If we knew $\theta$ and $\sigma^2$, then the prediction would be \begin{equation} p(y_{T+1}|x_{T+1},\theta,\sigma^2) = \textsf{N}(y_{T+1}|f(x_{T+1},\theta),\sigma^2) . \end{equation} The "prediction" is a probability distribution, from which we can compute the mean or the median or an interval --- whatever is appropriate for the problem.

Since we don't know $\theta$ and $\sigma^2$, we use the posterior distribution (which is what we do know) to integrate out our uncertainty: \begin{equation} p(y_{T+1}|y_{1:T},x_{1:T+1}) = \iint p(y_{T+1}|x_{T+1},\theta,\sigma^2)\,p(\theta,\sigma^2|y_{1:T},x_{1:T})\,d\theta\,d\sigma^2 . \end{equation} This "predictive distribution" is based on (1) the model, (2) what we have observed, and (3) the assumed covariate/condition for the "next" observation. It takes into account our uncertainty about the unknown model parameters and the unknown model error.