Obtain prediction interval for generalized linear model by simulation

When doing linear regression by LinearModel.fit, the "predict" command can return the prediction interval by setting the name-value pair 'Prediction' to 'observation'. However, such option doesn't exist when using the "predict" command on models obtained by GeneralizedLinearModel.fit. I read in few threads that this can be solved by simulation:

http://stats.stackexchange.com/questions/29962/method-to-estimate-the-prediction-interval-for-glm-and-negative-binomial-distrib

https://stat.ethz.ch/pipermail/r-help/2003-May/033165.html

But I am not sure about the exact procedure. Can anyone please help me on this?