Predictive distribution in the Frequentist setting Let $\{X_t\}_{t \in T}$ be a time series, such that $X \sim F(\theta)$, for some arbitrary distribution $F$. Based on the observed values $\{x_1, x_2, \cdots, x_{t-1}\}$, suppose that I want to predict $x_t$.
What should I do? My suggestion would be: estimate $\theta$, say, by $\hat{\theta}$ (e.g., MLE), $2)$ simulate sufficiently many points from $F(\hat{\theta})$, and $3)$ average values from step "$2)$" in order to obtain $\hat{x}_{t}$.
However, the above procedure has two problems (maybe more): $A)$ I am not accounting for the uncertainty in estimating $\theta$, and $B)$ I am obtaining a "deterministic point prediction", as opposed to a "probabilistic forecast" (in the form of a predictive distribution, which is what I want).
Assuming $X$ is absolutely continuous and, therefore, has a pdf $f_X$, I think I could address item $B)$ by arguing that $f_X(x|\hat{\theta})$ describes the predictive distribution. Is it correct? Even if it is true, I still have the problem described on item $A)$.
Also, under this same framework, how could I predict $\{x_{t+1}, x_{t+2}, \cdots, x_{t + p}\}$ based on the same original data set? In this case, the uncertainty must increase with $t$, right?

I think my question can be summarized by: "under the frequentist framework, how to make 'probabilistic forecasting' appropriately?".

Most forecasters silently ignore the uncertainty in parameter estimation, because (a) there are no nice closed form solutions as there are in OLS (using a $t$ distribution instead of a normal, in order to account for the variance of parameter estimates), and (b) to be honest, it often does not really matter all that much.
That said, you could of course try to get a handle on the estimation uncertainty of $\hat{\theta}$ (which, again, is in itself not easy to do in a time series context). Then get hopefully better predictive distributions by simulating many $\hat{\theta}_i$ according to the parameter distribution you assume, then for each $i$, drawing $x_i\sim F(\hat{\theta}_i)$. Then you can hope that the $(x_i)$ describe the true future density better than just a single $F(\hat{\theta})$.
This can then be iterated to get multi-step-ahead densities. Always assuming, of course, that you get a distributional handle on the corresponding parameter distributions. Any "reasonable" way of doing that should naturally yield larger variances for parameters that are farther out.
Note that usually, there is also uncertainty in just how $\hat{\theta}$ relates to previous observations - in an ARIMA setting, that would be uncertainty about the ARIMA orders. And similarly, assuming that we know the distributional family $F$ may also be a heroic assumption.
