First of all let me start by saying that I'm not used to using Python. Another thing is that you might want to think about your title again, more specific, the "correctly" term: There seems to be no evidence supporting that you can 100% accurately predict stock returns. Perhaps, you could look up the Efficient Market Hypothesis, more precise how the semi-strong form relates to GARCH models.
I guess you could use this method you're describing, however as far as I can see without actually calculating it myself, when you take the mean at the end you effectively kills the randomness of it. That is, when you generate X different point forecasts, they should in theory be distributed with the assumed distribution in the innovations (which are IID(0,1)), scaled with the (same) volatility forecast and mean, and thus, when taking the mean in the end, you simply obtain a number close to the forecasted conditional mean.
I would simply use a random generated value from the distribution you're assuming in the GARCH, for example a standard normal distribution. Remember that you assume $e \sim IID(0,1)$. You could look into a one-step ahead rolling forecast scheme and perhaps just check to see how your rolling forecast compare with real observed returns.
Another thing you could do is look up parametric density forecasting, in which you practically forecasts the volatility and then scale the assumed distribution of the innovations with this volatility.