scipy.stats.johnsonsu 中的 a 和 b 参数是什么？

Question

我正在尝试将Johnson SU分布拟合到一组经验丰富的标准普尔 500 指数回报中。我的理解（免责声明：不是数学专家）是这个分布包含了第三和第四时刻（偏斜和峰度）。除了loc（均值）和scale（标准差）johnsonsu之外，a还有两个额外的参数，和b。但是这些参数的顺序和规格令人困惑。

这就是我的困惑源于：如果我将回报纳入 SPDR S&P 500 ETF Trust (SPY)，我会得到以下经验统计数据：

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

现在，如果我对这个经验数据拟合正态分布，.fit应该返回loc和scale参数。（正态分布所需的一切。）检查：

import scipy.stats as scs

normmean, normstd = scs.norm.fit(r)
print(np.allclose(normmean, mean))
print(np.allclose(normstd, std))

True
True

但不太清楚返回的是什么scs.johnsonsu.fit：

print(scs.johnsonsu.fit(r))
(0.098009661042083682, 1.022060362199493, 0.0013471690867203458, 0.0072653444313926403)

这些应该是分布的四个参数：xi, gamma, delta, lam。

但我不能让他们回到经验平均值，这应该是：

IE：

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

和

mean, var, skew, kurt = scs.johnsonsu.stats(loc=xi, scale=lam, 
                                            a=gamma, b=delta, moments='msvk')

得到一堆NaNs。

Answer 1

它们是Johnson SU 的参数。请记住，您获得的样本均值与分布均值不同。这是平均值的表达式

这是方差的表达式：

在您的代码中，ξ 为loc，λ 为scale，γ 为a，δ 为b。sinh ^-1 (x) 等于 log(x + sqrt(1 + x ² ))。

因此检查拟合的返回值，为所有四个参数赋值，然后计算分布均值并与样本均值进行比较。如果有效，请重复练习以获得差异

更新

我试过你的代码并建议检查均值和方差，效果很好，请检查下面

import sys
import math

from pandas_datareader.data import DataReader as dr
import scipy.stats as scs

def read_data():
    return dr('SPY', 'google', start='2000')['Close'].pct_change().dropna()

def johnsonsu_mean(a, b, loc, scale):
    """
    Johnson SU mean according to https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution
    """
    v = loc - scale * math.exp(0.5 / b**2) * math.sinh(a/b)
    return v

def johnsonsu_var(a, b, loc, scale):
    """
    Johnson SU variance according to https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution
    """
    t = math.exp(1.0 / b**2)
    v = 0.5*scale**2 * (t - 1.0) * (t * math.cosh(2.0*a/b) + 1.0)
    return v

def johnsonsu_median(a, b, loc, scale):
    """
    Johnson SU median according to https://en.wikipedia.org/wiki/Johnson%27s_SU-distribution
    """
    v = loc + scale * math.sinh(-a/b)
    return v

def main(r):
    sample_mean, sample_med, sample_var, sample_std, sample_skew, sample_kurt = r.mean(), r.median(), r.var(0), r.std(0), r.skew(), r.kurt()

    a, b, loc, scale = scs.johnsonsu.fit(r) # fit the data and get distribution parameters back

    # distribution mean and variance according to SciPy
    dist_mean = scs.johnsonsu.mean(a, b, loc, scale)
    dist_med  = scs.johnsonsu.median(a, b, loc, scale)
    dist_var  = scs.johnsonsu.var(a, b, loc, scale)

    # distribution mean, var vs sample ones
    print("{0} {1}".format(sample_mean, dist_mean))
    print("{0} {1}".format(sample_med, dist_med))
    print("{0} {1}".format(sample_var, dist_var))
    print("")

    # distribution mean and variance according to Wiki vs SciPy
    print("{0} {1}".format(dist_mean, johnsonsu_mean(a, b, loc, scale)))
    print("{0} {1}".format(dist_var, johnsonsu_var(a, b, loc, scale)))
    print("{0} {1}".format(dist_med, johnsonsu_median(a, b, loc, scale)))

if __name__ == "__main__":
    r = read_data()
    main(r)

    sys.exit(0)

并产生输出：

0.00028012130615107805 0.00021391570000183283  
0.0005697194131890626 0.0006458197694718355    
0.00014415554662672425 0.00015479059187195545  

0.00021391570000183283 0.00021391570000541633  
0.00015479059187195545 0.00015479059186527505  
0.0006458197694718355 0.0006458197694718355