What is the Skewness of a Geometric Brownian Motion? Consider a GBM :
  $$S(t) = S(0)\exp\left({(\mu-\frac{1}{2}\sigma^2) t + \sigma W_t}\right)$$
  $$d\log S(t) = (\mu-\frac{1}{2}\sigma^2) t + \sigma dW_t$$
   $$\frac{d S(t)}{S(t)} = \mu t + \sigma dW_t$$
   where $W_t$ is the value of random Brownian Motion with difference distributed standard normal.
I have seen that 
$$\mathbb{E}(S_t)= S_0e^{\mu  t}$$
$$\operatorname{Var}(S_t)= S_0^2e^{2\mu t} \left( e^{\sigma^2 t}-1\right)$$
And I know that the median is $$S(0)e^{(\mu-\frac{1}{2}\sigma^2)t}$$, But I'm looking for an expression for the skewness (and kurtosis) if there is one.
I think the probability density function of a ''St'' is:
$f_{S_t}(s; \mu, \sigma, t) = \frac{1}{\sqrt{2 \pi}}\, \frac{1}{s \sigma \sqrt{t}}\, \exp \left( -\frac{ \left( \ln s - \ln S_0 - \left( \mu - \frac{1}{2} \sigma^2 \right) t \right)^2}{2\sigma^2 t} \right).$
Many thanks!
 A: To calculate the skewness (in fact, every order moments of GBM), we need Ito's lemma, https://en.wikipedia.org/wiki/It%C3%B4%27s_lemma, which states that for  an Ito diffusion
$$
dS_t =\mu_tdt+\sigma_tdW_t
$$
and a twice-differentiable scalar function $F(t, S)$, we have
$$
\begin{align}
F(t,S_t)=F(0, S_0)&+\int ^t _0 F_1 (s, S)ds\\
&+\int ^t _0 F_2 (s, S)\mu_sds\\
&+\int ^t _0 F_2 (s, S)\sigma_sdW_s\\
&+\dfrac {1}{2}\int ^t _0 F_{22} (s, S)\sigma_s ^2ds\\
\end{align}
$$
for all $t$.
These four integrals come from Taylor expansion and the property of Brownian motion that $(dW_t)^2=dt$.
For geometric Brownian motion, we see that
$$
\mu_t =\mu S_t, \sigma_t =\sigma S_t
$$
moreover, to calculate skewness, we let $F(t, S)=S^3$, which doesn't depend on $t$. Then Ito's lemma gives 
$$
S^3_t=S_0^3 + \int_0^t (3\mu +3\sigma^2)S_s^3ds+ \int_0^t 3\sigma S_s^3dW_s
$$
taking expectations, we get
$$
\mathbb E S^3_t= S_0^3+ \int_0^t (3\mu +3\sigma^2)\mathbb E S_s^3ds
$$
because the expectation of stochatic integral is always zero (provided $S_s^3$ is integrable).
Now let $h(t)=\mathbb E S^3_t$, then the above equation is
$$
h(t)=S_0^3+(3\mu +3\sigma^2)\int_0^t h(s)ds
$$
it's obvious:
$$
h(t)=S_0^3 e^{(3\mu +3\sigma^2)}=\mathbb E S^3_t
$$
Now we have all gradients for skewness: $\mathbb E S_t, \operatorname{Var} S_t, \mathbb E S^3_t$:
$$
\mathbb E [S_t-\mathbb E S_t]^3=\mathbb E S_t^3-3\mathbb E S_t\operatorname{Var} S_t-(\mathbb E S_t)^3
$$
and the result follows:
$$
S_0^3 e^{3\mu t}[e^{3\sigma^2 t}-3e^{\sigma^2 t}+2]
$$
