title: What are the a and b parameters in scipy.stats.johnsonsu?
body: <p>I am toying around with fitting a <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.johnsonsu.html" rel="nofollow noreferrer">Johnson SU</a> distribution to an empirical set of S&amp;P 500 Index returns.  My understanding (disclaimer: not a math guy) is that this distribution incorporates the third and fourth moments (skew &amp; kurtosis).  In addition to <code>loc</code> (the mean) and <code>scale</code> (the std deviation), <code>johnsonsu</code> takes two additional parameters, <code>a</code> and <code>b</code>.  But the order and specification of these parameters is confusing.</p>
  
  <p>Here's where my confusion stems from: if I pull in returns to SPDR S&amp;P 500 ETF Trust (SPY), I get the following empirical stats:</p>
  
  <pre><code>from pandas_datareader.data import DataReader as dr
  r = dr('SPY', 'google', start='2000')['Close'].pct_change().dropna()
  mean, var, std, skew, kurt = r.mean(), r.var(0), r.std(0), r.skew(), r.kurt() # ddof = 0
  # mean: 0.00027732907268771364
  # var: 0.00014416720067485022
  # std: 0.012006964673673785
  </code></pre>
  
  <p>Now if I fit a normal distribution to this empirical data, <code>.fit</code> is supposed to return the <code>loc</code> and <code>scale</code> parameters.  (All that is needed for a normal distribution.)  That checks out:</p>
  
  <pre><code>import scipy.stats as scs
  
  normmean, normstd = scs.norm.fit(r)
  print(np.allclose(normmean, mean))
  print(np.allclose(normstd, std))
  
  True
  True
  </code></pre>
  
  <p>But it is less clear what is returned by <code>scs.johnsonsu.fit</code>:</p>
  
  <pre><code>print(scs.johnsonsu.fit(r))
  (0.098009661042083682, 1.022060362199493, 0.0013471690867203458, 0.0072653444313926403)
  </code></pre>
  
  <p>These should be the <a href="https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution" rel="nofollow noreferrer">four parameters</a> to the distribution: <em>xi, gamma, delta, lam</em>.</p>
  
  <p>But I cannot get them to back into the empirical mean, which should be:
  <a href="https://i.stack.imgur.com/hGHDJ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/hGHDJ.png" alt="enter image description here"></a></p>
  
  <p>I.e.: </p>
  
  <pre><code>def johnsonmean(gamma, xi, delta, lam):
      mean = xi - lam * np.exp(delta ** -2 / 2) * np.sinh(gamma / delta)
      return mean
  gamma, xi, delta, lam = scs.johnsonsu.fit(r) # correct order?
  print(johnsonmean(gamma, xi, delta, lam))
  -inf
  </code></pre>
  
  <p>And</p>
  
  <pre><code>mean, var, skew, kurt = scs.johnsonsu.stats(loc=xi, scale=lam, 
                                              a=gamma, b=delta, moments='msvk')
  </code></pre>
  
  <p>gets a bunch of <code>NaN</code>s.</p>