Arbitrary prediction interval with random independent variables I have a equation with the parameters determined from multiple linear regression:
$$
Y = \beta_0 + X_1 \beta_1 + X_2\beta_2
$$
I would like to forecast the distribution of $Y$ numerically, in other words a histogram of the values $Y$ can take on in the future. From the histogram I can then deduce arbitrary prediction intervals. 
Both $X_1$ and $X_2$ are random variables with arbitrary distributions. For example $X_1$ could be uniform and $X_2$ could be normal. 
My approach
Use a Monte Carlo type approach. Randomly sample $X_1$, $X_2$, $\beta_1$ and $\beta_2$ and calculate $Y$, plot histogram. 


*

*Is this a reasonable approach? 

*If yes, then what is a reasonable distribution for $\beta_1$ and $\beta_2$?

*If no, how would I go about computing such a distribution. 


Many thanks in advance. 
Edit
Put another way: 
I would like to forecast $Y$ using forecast values of $X_1$ and $X_2$.  Instead of just producing a point estimate and 95% prediction interval I would like to produce a histogram which I can then use to get arbitrary quantiles of $Y$. I would like the forecast to take into account the fact that $X_1$ and $X_2$ have a given distribution. 
My simplest possible question is: 
How do I produce a histogram of forecast values of $Y$ that takes account of the fact that my regression is only an estimate and that my forecasts of $X_1$ etc. are random with a given distribution?
 A: What's missing from your approach, as @whuber notes in comments, is that your model really is;
$$
Y = \beta_0 + X_1 \beta_1 + X_2\beta_2 + \epsilon
$$
where $\epsilon$ is the residual error.
So you have to consider not only the distributions of the $\beta_i$ but also of $\epsilon$.
Your regression analysis probably has enough information to make some reasonable estimates. Depending on what you are trying to accomplish, you might just get by with using the standard errors of the $\beta_i$ for the distributions of coefficients and the residual mean-square error for the distribution of $\epsilon$, taking them as normally distributed. You are probably better off exploring the distribution of residual errors as functions of $X_1$ and $X_2$ jointly to see if such an assumption might get you into trouble, as @whuber indicated.
Instead of assuming distributions of the $\beta_i$, you could also consider repeating your analyses on bootstrapped samples of your original data, to determine their quasi-empirical distributions.
Provided that you specify the nature of your forecasts, I don't see that there should be an issue with reporting such histograms.
