Why does the sim function in Gelman's arm package simulate sigma from inverse chi square? In getMethod(arm::sim, "lm"), the source code shows that $\sigma$ is simulated from inverse chi square:
for (s in 1:n.sims) {
            sigma[s] <- sigma.hat * sqrt((n - k)/rchisq(1, n - 
                k))
            beta[s, ] <- MASS::mvrnorm(1, beta.hat, V.beta * 
                sigma[s]^2)
        }

This makes sense if one is doing Bayesian linear regression with the default uniform prior for $(\beta, \log \sigma)$. But if I'm running frequentist MLE regression, isn't it wrong to use arm::sim to do simulation? Indeed, for MLE, the asymptotic variance of $\sigma^2$ is normal instead.
 A: Using a uniform prior is equivalent to the MLE estimate. When you say that the asymptotic variance of $σ^{2}$ is normal, that's not quite correct. The quantity $(n-1)×s^{2}/σ^{2}$ is distributed as $χ^{2}_{n-1}$, with $s$ the sample standard deviation, $σ$ the true population standard deviation, and $n$ the sample size. With large $n$, a $χ^{2}$ distribution is essentially normal, which might be the source of the confusion here. Extending this for the case of residual variance in a linear model, the quantity $(n-k)×σ̂^{2}/σ^{2}$ is distributed as $χ^{2}_{n-k}$, with $k$ the number of predictors, $σ̂$ the estimated residual standard deviation from the model, and other terms defined as above.
So, the arm::sim is drawing sigma values from a standard chi square distribution (note, not an inverse chi square distribution as the title states), then scaling the values up by the estimated $σ̂$ value so that the distribution is on the same scale as the sample data. This is the the appropriate frequentist approach to simulating samples from a linear model.
