Is there a way to generate non-normal correlated random variates?

This is my first post so excuse me for any inconveniences.

My goal is to simulate asset return timeseries, which typically are fat-tailed. Therefore, I must look further than generating correlated normal variates. I came across this package that basically does exactly what I want. However, it is for python 2.x.x and I'm on python 3.8.

from connorav import CorrelatedNonNormalRandomVariates
returns = pd.read_csv(r'datasets/eurostoxx_share_returns.csv', sep=';', parse_dates=True, index_col='Date')
nsim = 10
corr = returns.corr()
corr = posdef.nearestPD(returns.cov())
stats = pd.DataFrame(data={'mean': returns.mean(),
                           'std': returns.std(),
                           'skew': returns.skew(),
                           'kurt': returns.kurt()
                           }).to_numpy()
corr = corr.to_numpy()
rv = CorrelatedNonNormalRandomVariates(stats,corr,nsim)

I received an error because python 2 would return a list upon map() where python 3 does not. So I adjusted in correl_rv.py:

self.moments = moments
self.correlations = correlations
self.distributions = map(MSSKDistribution,moments.tolist())

to

self.moments = moments
self.correlations = correlations
self.distributions = list(map(MSSKDistribution,moments.tolist()))

Yet now I receive the following output:

~\Anaconda3\envs\python\lib\site-packages\connorav\correl_rv.py in __init__(self, moments, correlations, num_samples, method)
     21         self.moments = moments
     22         self.correlations = correlations
---> 23         self.distributions = list(map(MSSKDistribution,moments.tolist()))
     24 
     25         self.generate(num_samples,method)

~\Anaconda3\envs\python\lib\site-packages\connorav\distribution.py in __init__(self, mean, std, skew, kurt)
     12         self.skew = skew
     13         self.kurt = kurt
---> 14         self.fit()
     15 
     16     def fit(self):

~\Anaconda3\envs\python\lib\site-packages\connorav\distribution.py in fit(self)
     16     def fit(self):
     17 
---> 18         if abs(self.skew) < NORMAL_CUTOFF and abs(self.kurt) < NORMAL_CUTOFF:
     19             # It is hard to solve the johnson su curve when it is very close
     20             # to normality, so just use a normal curve instead.

TypeError: bad operand type for abs(): 'NoneType'

The data I use are daily stock returns from all 50 EuroSTOXX firms over 2000-2021.

Thanks to you all!