We develop an approach to approximate the multivariate distribution of time-aggregated stock returns in the generalized autoregressive conditional heteroscedasticity (GARCH) context. The approach yields a single-time-step simulation procedure as opposed to the multiple-time-step simulation currently required in this context. For this purpose, the exact moment formulas for the time-aggregated return under a quadratic GARCH process are combined with multivariate nonnormal simulation procedures using as inputs the first four moments and correlation structure of the unknown target distribution. Estimation and simulation results are presented for a portfolio of thirty stocks from the Dow Jones Industrial Average index. The results reveal that the proposed simulation method can generate random numbers with moments and correlations agreeing with the targets. Using value-at-risk computations for different horizons and probabilities, we show that the percentiles of portfolios' return distributions computed with the proposed approach provide good approximations to benchmark values obtained from a multistep simulation.