Correct way to apply model uncertainty in bootstrapping Monte Carlo I have a simple linear model that predicts sales as a function of a temperature-related variable. 
For the sake of simplicity, let’s assume that the model is 
$$\hat{Y} = a+b \cdot T$$
I train this on a set of historical values $\{T_i, Y_i\}$. The estimated coefficients $a$ and $b$ have their own statistical uncertainty 
and the residuals $e_i = (\hat{Y}_i - Y_{i})$ are gaussianly distributed with $N(0, \sigma^2_e)$ with
$$\sigma^2_e = \sum{(\hat{Y}-Y_i})^2/(N-2)$$   (the standard error of the estimate)
Notice that this variance is intrinsic in the data, in the sense that for the same temperature I can have a range of $Y$ values. 
This can be imputed to missing explanatory variables in my model or genuinely random behavior.
Then I want to forecast the number of sales, taking into account the natural variability of weather conditions.   
Because I don’t have a model that accurately describes the p.d.f. of the temperature, I bootstrap this using a long series of temperature data.
I simulate thousands of synthetic weather days in this way and for each one of them I calculate an estimated number of sales with my formula. 
On top of the weather variability, there is my model uncertainty. I want that my simulations correctly capture the model uncertainty mentioned above: for the same temperature, I should get different calculations for $Y$. I smear gaussianly the predictions.
Question: The question is what is the correct sigma to use when I apply the Gaussian smearing.
I’m considering three options but I don’t know which one is correct.
1. Use the standard error of the estimate, $\sigma_e$
2. Compute the uncertainty of the prediction, $\sigma_{pred}$ applying the error propagation rules 
to the stat uncertainties of the coefficients and the input temperatures used in the simulation
3. Sum in quadrature $\sigma_e$ and $\sigma_{pred}$
What is the correct thing to do in this case?
 A: When you say bootstrap Monte Carlo I assume you really meant bootstrap.  Monte Carlo simulation is for when you know the true underlying (population) model.
You did not state what types of model uncertainties pertain.  You wrote more about parameter estimate uncertainty and didn't mention entertaining any other distributions, transformations, or predictors.  Assuming you only need to deal with parameter (regression coefficient) uncertainty, all the mechanisms of standard regression modeling apply and your solution is easy with a little study of statistics.  Get the prediction interval and a confidence interval (e.g. 0.95 confidence interval) for the prediction.  Note this is a much wider interval for the usual one which is for estimating long-term averages ($E(Y|X)$).  Confidence limits for individual $Y$ are dominated by the estimated residual variance of $Y|X$ but also contain a variance component describing the difficulty of estimating the regression parameters.  It is only this component to converges to zero as $n \rightarrow \infty$.
Your comment about imputing missings is very hard to understand.
