This can be done using the sinh-arcsinh transformation from

> Jones, M. C. and Pewsey A. (2009). [Sinh-arcsinh distributions](http://biomet.oxfordjournals.org/content/96/4/761.abstract). Biometrika 96: 761–780.

The transformation is defined as 

$$H(x;\epsilon,\delta)=\sinh[\delta\sinh^{-1}(x)-\epsilon], \tag{$\star$}$$

where $\epsilon \in{\mathbb R}$ and $\delta \in {\mathbb R}_+$. When this transformation is applied to the normal CDF $S(x;\epsilon,\delta)=\Phi[H(x;\epsilon,\delta)]$, it produces a distribution where the parameters $(\epsilon,\delta)$ control skewness and kurtosis, respectively (Jones and Pewsey, 2009), in the sense of [van Zwet (1969)](http://www.amazon.co.uk/Convex-transformations-random-variables-Mathematical/dp/B0007J048W). In addition, if $\epsilon=0$ and $\delta=1$, we obtain the original normal distribution. See the following R code.

    fs = function(x,epsilon,delta) dnorm(sinh(delta*asinh(x)-epsilon))*delta*cosh(delta*asinh(x)-epsilon)/sqrt(1+x^2)
    
    vec = seq(-15,15,0.001)
    
    plot(vec,fs(vec,0,1),type="l")
    points(vec,fs(vec,1,1),type="l",col="red")
    points(vec,fs(vec,2,1),type="l",col="blue")
    points(vec,fs(vec,-1,1),type="l",col="red")
    points(vec,fs(vec,-2,1),type="l",col="blue")
    
    vec = seq(-5,5,0.001)
    
    plot(vec,fs(vec,0,0.5),type="l",ylim=c(0,1))
    points(vec,fs(vec,0,0.75),type="l",col="red")
    points(vec,fs(vec,0,1),type="l",col="blue")
    points(vec,fs(vec,0,1.25),type="l",col="red")
    points(vec,fs(vec,0,1.5),type="l",col="blue")

Therefore, by choosing an appropriate sequence of parameters $(\epsilon_n,\delta_n)$, you can generate a sequence of distributions/transformations with different levels of skewness and kurtosis and make them look as similar or different to the normal distribution as you want.

The following plot shows the outcome produced by the R code. For *(i)* $\epsilon=(-2,-1,0,1,2)$ and $\delta=1$, and *(ii)* $\epsilon=0$ and $\delta=(0.5,0.75,1,1.25,1.5)$.

![enter image description here][1]

![enter image description here][2]

Simulation of this distribution is straightforward given that you just have to transform a normal sample using the inverse of $(\star)$.

$$H^{-1}(x;\epsilon,\delta)=\sinh^{-1}[\delta^{-1}\sinh(x)+\epsilon]$$


  [1]: https://i.sstatic.net/cbQQR.jpg
  [2]: https://i.sstatic.net/d1SWf.jpg