Correct procedure to simulate prediction intervals

Let $T_i$ be the survival time of patient $i$ and $x_i$ (a single covariate). I assume a log-normal accelerated failure time model: $$ y_i = \log T_i = \beta_0 + \beta_1 x_i + \sigma \varepsilon_i, \qquad \varepsilon_i \overset{\text{iid}}{\sim} \mathcal N(0,1), \quad i = 1, \dots, N, $$ $$y_i \mid x_i \sim \mathcal N(\beta_0 + \beta_1 x_i,\ \sigma^2)$$ I want to know how to construct a prediction interval for the (log) survival time of a single new patient, and for the average (log) survival time of a group of $n$ new patients who all have the same covariate value $x_0$ . My understanding is that for a new patient with covariate value $x_0$ , the actual outcome is $y_0 = \beta_0 + \beta_1 x_0 + \sigma\varepsilon_0$ , so $E[y_0 \mid x_0] = \beta_0 + \beta_1 x_0$ . The mean prediction $\hat y_0 = \hat\beta_0 + \hat\beta_1 x_0$ is an estimate of $E[y_0 \mid x_0]$ . This is my understanding on the procedure for prediction intervals for both individual and groups: Prediction interval for the log survival time of a single new patient $i$ with covariate $x_0$ : Step 1. Fit model by MLE to training data to get $\hat\theta = (\hat\beta_0, \hat\beta_1, \log\hat\sigma)$ and the estimated covariance matrix $\widehat{\operatorname{Cov}}(\hat\theta)$ . Step 2. Choose the number of simulations $S$ (e.g. $S = 10{,}000$ ). Step 3. For simulation $s = 1$ , draw parameters from the asymptotic distribution of the MLE: $$ \theta^{(s=1)} = \big(\beta_0^{(s=1)}, \beta_1^{(s=1)}, \log\sigma^{(s=1)…