I am going to supply an answer that is quite similar to SRKX's (which is very very good) because I want to discuss in more detail a few important things. First, you cannot use a stochastic volatility model for the SDE that you've provided as that's GBM with constant diffusion. However, based on what you've said it's obvious you wish to model a discretized version of the following process:
$$\boxed{dS(t) = S(t)(\mu dt + \sigma(t)dW(t))}$$
where $W(t)$ is a standard scalar Brownian motion under the real-world probability measure, $\mu$ is its constant drift and $\sigma(t)$ is its volatility process which we assume to follow GARCH dynamics.
Doing what we always do most of the time in finance research, I'm going to define returns as $R(t):= \ln{\frac{S(t)}{S(t-1)}}$. In the univariate setting, the mean equation of a ARCH/GARCH type model generally follows:
$$\boxed{R(t) = \mu + \sum_{i=1}^N a_i R(t-i) + b \sigma(t) + \mathbf{c}\mathbf{X(t)} + \sigma(t) (W(t)-W(t-1))}$$
where $a_i$ are the AR(i) coefficients, $b$ is the coefficient of the garch-in-mean, $\mathbf{c}$ is a row vector of coefficients to some exogenous variables represented by column vector $\mathbf{X(t)}$. However, given the process that you've specified, we are restricted to a very specific mean equation:
$$\boxed{R(t) = \mu + \sigma(t)(W(t)-W(t-1)) }$$
where $W(t)-W(t-1) \overset{d}{=} W(1) \sim N(0,1)$. This restriction is not automatically troubling as it is prevalent when modelling DCC-fGARCH (time-varying correlations with family GARCH). However most studies that look to things such as bivariate variance spillovers (e.g., contagion research) or that are doing a study where univariate GARCH is involed (e.g., exchange rate exposure research where market returns and exchange rate returns are exogenous variables in the mean equation) will find this process to be inappropriate for their usage.
Ignoring the empirical issues with the mean equation that follow from the proposed SDE, we can see that the mean equation representation follows from the fact that, at $\mathcal{F}_{0}$ and under real-world $\mathbb{P}$:
$S(t) = e^{\mu t + \sigma(t) W(t)}$
$\therefore \ln{\frac{S(t)}{S(t-1)}} = \mu + \sigma(t) (W(t) - W(t-1))$
$\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \overset{d}{=}\mu + \sigma(t) W(1) $
Luckily, the variance equation is unconstrained, and we can use the GARCH model whose process is defined here. For a discretized econometric representation of GARCH(1,1) we have, as SRKX lays out;
$$\boxed{\sigma_t^2 = \omega + \alpha (W(t-1)-W(t-2))^2\sigma(t-1)^2 + \beta \sigma(t-1)^2}$$
So you should then proceed with your plan of discretization while using the correct mean equation that I boxed. This is the default in R packages ccgarch
, rugarch
and fGarch
, so it's your lucky day!