How to simulate non-gaussian stochastic paths (Edited to be clearer)
I am trying to replicate simulating Geometric Brownian Motion (GBM) but instead of the stochastic increment following a normal distribution, I would like it to follow a distribution of my choice.
With standard GBM, if we have a standard deviation of 10% for a full period then we have the following results if we were to simulate a load of paths:

*

*the full period distribution of those paths will be normal with a std of 10%

*if we split the period up into 100 increments, any given increment will exhibit a normal distribution with a std of 1%

*in general for a given increment with std of 1%, a period of n increments will be normally distributed with a std of 1% * sqrt(n)

In other words, we simulate each increment from a normal distribution with a std of x% and as the paths continue through time they continue to display a normal distribution but with the standard deviation scaling with the square root of time i.e. x% * sqrt(n)
I would like to achieve the same scaling but with a different (fatter tailed distribution) - so the variance scales and the distribution remains the distribution I want. The problem I am encountering is that simply changing the distribution of the increment to my new distribution is not good enough. This is because if I generate 100 increments from a new distribution then the sum of these increments becomes normal - and so the distribution of my paths heads toward normal very quickly due to CLT.
In code, standard GBM can be written as the following in python:
import numpy as np

def generate_paths(spot, drift, sigma, years, points, sims):

    # create numpy PCG PRNG
    rng = np.random.default_rng(seed=123)
    # create random numbers
    dW = rng.normal(size=(sims, int(points * years)))
    # get our scaled random deviations
    scaled_dW = dW * (sigma / points**0.5)
    # add on the drift
    daily_devs = (drift / points) + scaled_dW
    # cumsum them as log returns additive
    cum_rets = np.cumsum(daily_devs, axis=1)
    # add in 0 at start for starting price
    cum_rets = np.insert(cum_rets, 0, 0, axis=1)
    # create price series
    pxs = spot * np.exp(cum_rets)
    # return the goodies
    return daily_devs, pxs

Instead of this I'd like to change:
dW = rng.normal(size=(sims, int(points * years)))

for something like:
dW = rng.NEW_DIST(size=(sims, int(points * years)))

and this new distribution to persist s.t. CLT doesn't mean that the distribution of my paths heads towards normal quickly.
Here's an example of simulations I have tried to run - left to right shows the distribution of the paths after 1, 10, 100 and 250 increments where 250 is the max. As you can see after 1 increment we have the distribution I want - fat tailed. As soon as we start to aggregate and move through time we start to recover the normal distribution.

My initial thoughts are along the lines of:

*

*any sum of independent random increments with finite variance will converge to the normal distribution due to CLT

*perhaps this means the only way to circumvent is through introducing some correlation into the individual increments

*another solution might be to choose individual increments that follow an infinite variance process but this might be at odds with the data

 A: Posting as an answer as too long for a comment:
The reason you're seeing the central limit theorem crop up here is because your returns at each time point are independent.
I think what you want to do and what you're saying you want to do aren't quite the same.
If you simply had a time series of distributions that you wanted to match, then you could use the method Dupire used in order to determine the local volatility surface. This allows you do create a function \hat{\sigma}(S_t,t) such that the diffusion process:
$$ \frac{\mathrm{d}S}{S} = \mu(t) \mathrm{d}t + \hat{\sigma}(S_t,t)\mathrm{d}W_t$$
produces a probability distribution as a function of time that matches your requirement.
The above is simply the Local Volatility model, which sits on top of GBM and has $\mathrm{d}W_t$ distributed normally. What you're looking to do is essentially the same, but with an added twist - you're augmenting the returns. This will meaningfully change the maths, and mean much of the derivations won't apply, but that doesn't really matter, as you can just back solve. What you're essentially looking to do is the following:
$$ \frac{\mathrm{d}S}{S} = \mu(t) \mathrm{d}t + \hat{\sigma}(S_t,t|\bar{x})f(\mathrm{d}W_t)$$
or
$$ \frac{\mathrm{d}S}{S} = \mu(t) \mathrm{d}t + \hat{\sigma}(S_t,t|\bar{x})\mathrm{d}\tilde{W}_t$$
where $f(\mathrm{d}{W}_t) = \mathrm{d}\tilde{W}_t \nsim \mathcal{N}(0,1)$, where $\hat{\sigma}(S_t,t|\bar{x})$ is your local volatility for your given formulation of $f$, and $\bar{x}$ are you volatility model parameters.
IMO, as soon as you start moving into this realm, then the maths gets very difficult, and it's easier to instead just write your montecarlo code and then throw the whole thing through a solver to obtain some fit params for your local volatility formulation.
