How to find the best fitting parametric distribution for an empirical dataset (stock returns)?
First of all, if you want to find the best distribution that fits your data you just iteratively fit your data to the longlist of distributions. Scipy supports most of them. After fitting, you can either use KS-test to find which distribution fitted best or you can use fit error to decide. This solution does what you want, also other solutions in that post are very good approaches to your problem.
Also, if you are sure about your distribution's being Johnson-SU distribution, to find the parameters, use Scipy's fit function, which will return you the parameters that best represent your data.
a, b, loc, scale = scipy.stats.johnsonsu.fit(data)