Please let me know where I have been mistaken!
Let the SDE satisfied by the GBM $S(t)$ be $$ \frac{dS(t)}{S(t)} = \mu dt + \sigma dW(t). $$
Then, the underlying BM $X(t)$ will satisfy $$ dX(t) = \left( \mu - \frac{1}{2} \sigma^2 \right)dt + \sigma dW(t). $$
To simulate the GBM for times $t_0 < t_1 < \ldots < t_n$, generate $n$ iid $\mathcal{N}(0,1)$ RVs, $Z_i$, $\quad i = 1,2,\ldots, n$ and set
$$ S(t_{i+1}) = S(t_i)\exp\left(\left( \mu - \frac{1}{2} \sigma^2 \right)(t_{i+1} -t_i) + \sigma \sqrt{t_{i+1} - t_i} Z_{i+1} \right). $$ Then $$ \frac{S(t)}{S(0)} = \exp\left(\underbrace{\left(\mu - \frac{1}{2} \sigma^2\right) t + \sigma \sqrt{t}Z}_{\mathcal{N}\left(\left(\mu - \frac{\sigma^2}{2}\right)t, \sigma^2 t\right)} \right) $$ and so $$ \log \frac{S(t)}{S(0)} \sim \mathcal{N}\left(\left(\mu - \frac{\sigma^2}{2}\right)t, \sigma^2 t\right). $$
I think I'm missing something because when I use this mean and variance (in the equation directly above) for testing the normality of the log returns ($\log \frac{S(t)}{S(0)}$), I get ridiculous answers.
To be specific, with $S_0 = 20$, $\mu = 2$, $\sigma^2 = 1$ and partitioning $[0,1]$ into 100 subintervals, generating the GBM at these 100 points gives a range of values from 15.399 to 97.1384 for the $S(t_i)$. Then the log returns, $\log S(t_i)/S_0$, range from -0.26143 to 1.5804. The means I'm using for each $\log S(t_i)/S_0$ are (in increasing order of $i$) $0.015, 0.03, 0.045, \ldots, 1.5$ and the variances are $0.01,0.02, \ldots, 1$ since these are the parameters in the supposed normal distribution as described above. Finally, when these log returns are normalized using this mean and variance their values range from -2.7643 to 0.080404, which is clearly not $\mathcal{N}(0,1)$-distributed.
Thanks in advance!