Podemos mostrar que los momentos del Árbol Binomial coinciden con los momentos del modelo continuo para el caso en que elegimos un valor de probabilidad simétrico $p=0.5$.
Cambiaré ligeramente la notación (verás más adelante por qué: deliberadamente uso $\eta$ en lugar de $\mu$):
$$u=e^{\eta\Delta t+\sigma\sqrt{\Delta t}}, \qquad d=e^{\eta\Delta t-\sigma\sqrt{\Delta t}}$$
También, introduzcamos la siguiente notación: la madurez en la que estamos interesados será denotada como $T$ (en lugar de $t$), y $\Delta t$ puede expresarse como $\frac{T}{n}$, donde $n$ es el número de pasos en el árbol.
Por lo tanto, la dinámica para $S_t$ se da como:
$$S_T=S_0u^jd^{n-j}$$
Tomar el logaritmo hace las cosas mucho más fáciles:
$$\ln{\left(\frac{S_T}{S_0}\right)}=j\ln(u)+(n-j)\ln(d)=\\=j\left(\ln(u)-\ln(d)\right)+n\ln(d)=\\=2j\sigma\sqrt{\Delta t}+n\eta\Delta t-n\sigma\sqrt{\Delta t}=\\=j\left(2\sigma\sqrt{\frac{T}{n}}\right)+n\left(\eta\frac{T}{n}-\sqrt{\frac{T}{n}}\sigma\right)$$
Sabemos que $j\sim Bin(n,p)$, así que los momentos son:
\begin{align*} \tag{1} \mathbb{E}\left[\ln{\left(\frac{S_T}{S_0}\right)}\right]=\left(2\sigma\sqrt{\frac{T}{n}}\right)\mathbb{E}\left[j\right]+n\left(\eta\frac{T}{n}-\sqrt{\frac{T}{n}}\sigma\right)=\\=\left(2\sigma\sqrt{\frac{T}{n}}\right)(np)+n\left(\eta\frac{T}{n}-\sqrt{\frac{T}{n}}\sigma\right)=\\=2p(\sigma\sqrt{T}\sqrt{n})-\sigma\sqrt{T}\sqrt{n}+\eta T=\\=(2p-1)(\sigma\sqrt{T}\sqrt{n})+\eta T \end{align*}
La varianza es: \begin{align*} \tag{2} V\left(\ln{\left(\frac{S_T}{S_0}\right)}\right)=V\left(j2\sigma\sqrt{\frac{T}{n}}\right)=\\=V(j)\left(4\sigma^2\frac{T}{n}\right)\\=np(1-p)\left(4\sigma^2\frac{T}{n}\right)=\\=p(1-p)4T\sigma^2 \end{align*}
Si $p=0.5$, entonces los momentos anteriores son:
$$\mathbb{E}\left[\ln{\left(\frac{S_T}{S_0}\right)}\right]=\eta T \qquad V\left(\ln{\left(\frac{S_T}{S_0}\right)}\right)=T\sigma^2$$
Por lo tanto, si elegimos $\eta:=\mu-0.5\sigma^2$, entonces los momentos del modelo del Árbol Binomial coinciden con el modelo continuo siempre que $p=0.5$.
(la belleza es que si elegimos $p=0.5$ los momentos coinciden con el modelo continuo para cualquier $n$).
Entonces, con el escenario anterior, para un número (pequeño) finito de pasos "$n$" en el modelo del Árbol Binomial, la variable aleatoria $\ln{\left(\frac{S_T}{S_0}\right)}$ seguirá una distribución Binomial, pero su media y varianza coincidirán con el modelo GBM continuo.
¿Qué pasa con la convergencia del modelo de Binomial al modelo GBM para un $n$ grande? Aquí podemos mostrar convergencia en distribución:
Usando el TCL, podemos decir que para un $n$ grande:
$$j\xrightarrow{d}N(np,\sqrt{np(1-p)})$$
Entonces (usando los resultados (1) y (2)) podemos decir que para un $n$ grande:
\begin{equation} \boxed{\left( \frac{S_T}{S_0} \right)\xrightarrow{d}N\left(\sqrt{Tn}\sigma(2p-1)+\eta T,p(1-p)4T\sigma^2\right)} \end{equation}
Nuevamente, esto coincidirá con el modelo continuo siempre que $p=0.5$ y $\eta=\mu-0.5\sigma^2$.
Resumen:
- Para cualquier número finito ("pequeño") de pasos $n$ en el modelo del Árbol Binomial, los momentos coincidirán con el modelo continuo siempre que $p=0.5$ y $\eta = \mu-0.5\sigma^2$ (y bajo esta configuración, la variable $\ln\left(\frac{S_T}{S_0}\right)$ se distribuirá Binomialmente con estos momentos)
- Si, además, $n$ se vuelve grande, la distribución de la variable $\ln\left(\frac{S_T}{S_0}\right)$ convergerá de Binomial a Normal (con los mismos momentos)
Nota final: anteriormente observamos el caso donde $\mu$ es la deriva histórica. Esto nos da la libertad de elegir $p=0.5$. Si en cambio consideramos el modelo de riesgo neutro con $\eta=r-0.5\sigma^2$, entonces el parámetro $p$ está dado por $p:=\frac{e^{r\frac{T}{n}-d}}{u-d}$ y es mucho más difícil mostrar la convergencia.
Replicando tus resultados:
Trabajando línea por línea obtengo:
$$S_t=S_0e^{j\left(2\sigma\sqrt{\Delta t}\right)+n\left(\mu\Delta t-\sqrt{\Delta t}\sigma\right)}=\\=S_0\exp{\left(\mu t+2\left(j\sigma \sqrt{\Delta t} \right) - n\sigma\sqrt{\Delta t}\right)}=\\=S_0\exp{\left(\mu t+\sigma \sqrt{\Delta t}\left(2j - n \right) \right)}$$
Sustituyendo por $j=np+\sqrt{np(1-p)}z$, obtenemos:
$$S_0\exp{\left(\mu t+\sigma \sqrt{\Delta t}\left(2np+2\sqrt{np(1-p)}z - n \right) \right)}$$
Ahora evaluando los términos, obtenemos:
$$np=n(0.5+0.5\frac{\mu}{\sigma}\sqrt{\Delta t})=0.5n+0.5\frac{\mu}{\sigma}\frac{t}{\sqrt{\Delta t}}$$
Ahora el otro término:
$2\sqrt{np(1-p)}z=2\left(np-np^2\right)^{\frac{1}{2}}z=\\=2\left(0.5n+0.5\frac{\mu}{\sigma}\frac{t}{\sqrt{\Delta t}}-n(0.5 -0.5\frac{\mu}{\sigma}\sqrt{\Delta t})^2\right)^{\frac{1}{2}}z=\\=2\left(0.5n+0.5\frac{\mu}{\sigma}\frac{t}{\sqrt{\Delta t}}-n(0.25-0.5\frac{\mu}{\sigma}\sqrt{\Delta t}+0.25\frac{\mu^2}{\sigma^2}\Delta t)\right)^{\frac{1}{2}}z=\\=2\left(0.5n+0.5\frac{\mu}{\sigma}\frac{t}{\sqrt{\Delta t}}-0.25n-0.5\frac{\mu}{\sigma}\frac{t}{\sqrt{\Delta t}}+0.25\frac{\mu^2}{\sigma^2}t\right)^{\frac{1}{2}}z=\\=2\left(0.25n+0.25\frac{\mu^2}{\sigma^2}t\right)^{\frac{1}{2}}z$
Reemplazándolo todo (usando $\sigma \sqrt{\Delta t}*2np=\mu t)$):
$$S_t=S_0\exp{\left(2\mu t+2\sigma \sqrt{\Delta t}\left(0.25n+0.25\frac{\mu^2}{\sigma^2}t\right)^{\frac{1}{2}}z\right)}=\\=S_0\exp{\left(2\mu t+2\sigma \sqrt{\Delta t} \left(0.25n+0.25\frac{\mu^2}{\sigma^2}n \Delta t\right)^{\frac{1}{2}}z\right)}=\\=S_0\exp{\left(2\mu t+\sigma \sqrt{t} \left(1+\frac{\mu^2}{\sigma^2} \Delta t\right)^{\frac{1}{2}}z\right)}$$
Para replicar tu resultado completamente, debemos mostrar que:
$$\lim_{n\to\infty}\left(1+\frac{\mu^2}{\sigma^2} \Delta t\right)^{\frac{1}{2}}=0$$