python - 理解Python中赫斯特指数的广义公式

Question

如果我的问题不是很简洁，我深表歉意。我试图理解 QuantStart 中发布的估计赫斯特指数的通用公式:

https://www.quantstart.com/articles/Basics-of-Statistical-Mean-Reversion-Testing

这是代码:

from numpy import cumsum, log, polyfit, sqrt, std, subtract
from numpy.random import randn
def hurst(ts):
    lags = range(10,80)
    tau = [sqrt(std(subtract(ts[lag:], ts[:-lag]))) for lag in lags]
    poly = polyfit(log(lags), log(tau), 1)
    return poly[0]*2.0
print ('> Hurst Exponent:', '%.6s' % hurst(df['Close']))

我已经应用了这段代码并且它有效。但我不太明白这个公式到底是如何运作的。我有疑问，尤其是滞后。如果有人能给我一些信息来理解公式的这一部分，我将非常感激。

例如，如果我将滞后范围设置为 10 和 80，这是否意味着将在时间序列中的 10 和 80 个观测值之间计算自相关？

另一个重要问题是，我需要多少个观测值才能设置滞后范围？我应该至少有 80 个观察结果吗？

正如我之前提到的，如果有人能好心地阐明这一点，我将非常感激。我将非常感激。

谢谢!

Answer 1

Q : how many observations do I need to set up that lags range?

考虑到需要 200 个柱的滞后，观测值少于 200 个柱确实没有意义。

Q : if I set up the lags range in 10 and 80, does that mean that the autocorrelation will be calculated between the 10 and 80 observations in a time series?

绝对是的。这就是赫斯特的工作方式。

查看函数实现，取自 QuantFX 模块，按原样 v.4.13(Py2.7 在大多数地方不会造成麻烦，但任何 xrange() 都应该在 Py3.x 中替换为 range()):

def HurstEXP( ts = [ None, ] ):                                         # TESTED: HurstEXP()                Hurst exponent ( Browninan Motion & other observations measure ) 100+ BARs back(!)
            """                                                         __doc__
            USAGE:
                        HurstEXP( ts = [ None, ] )

                        Returns the Hurst Exponent of the time series vector ts[]

            PARAMETERS:
                        ts[,]   a time-series, with 100+ elements
                                ( or [ None, ] that produces a demo run )

            RETURNS:
                        float - a Hurst Exponent approximation,
                                as a real value
                                or
                                an explanatory string on an empty call
            THROWS:
                        n/a
            EXAMPLE:
                        >>> HurstEXP()                                        # actual numbers will vary, as per np.random.randn() generator used
                        HurstEXP( Geometric Browian Motion ):    0.49447454
                        HurstEXP(    Mean-Reverting Series ):   -0.00016013
                        HurstEXP(          Trending Series ):    0.95748937
                        'SYNTH series demo ( on HurstEXP( ts == [ None, ] ) ) # actual numbers vary, as per np.random.randn() generator'

                        >>> HurstEXP( rolling_window( aDSEG[:,idxC], 100 ) )
            REF.s:
                        >>> www.quantstart.com/articles/Basics-of-Statistical-Mean-Reversion-Testing
            """
            #---------------------------------------------------------------------------------------------------------------------------<self-reflective>
            if ( ts[0] == None ):                                       # DEMO: Create a SYNTH Geometric Brownian Motion, Mean-Reverting and Trending Series:

                 gbm = np.log( 1000 + np.cumsum(     np.random.randn( 100000 ) ) )  # a Geometric Brownian Motion[log(1000 + rand), log(1000 + rand + rand ), log(1000 + rand + rand + rand ),... log(  1000 + rand + ... )]
                 mr  = np.log( 1000 +                np.random.randn( 100000 )   )  # a Mean-Reverting Series    [log(1000 + rand), log(1000 + rand        ), log(1000 + rand               ),... log(  1000 + rand       )]
                 tr  = np.log( 1000 + np.cumsum( 1 + np.random.randn( 100000 ) ) )  # a Trending Series          [log(1001 + rand), log(1002 + rand + rand ), log(1003 + rand + rand + rand ),... log(101000 + rand + ... )]

                                                                        # Output the Hurst Exponent for each of the above SYNTH series
                 print ( "HurstEXP( Geometric Browian Motion ):   {0: > 12.8f}".format( HurstEXP( gbm ) ) )
                 print ( "HurstEXP(    Mean-Reverting Series ):   {0: > 12.8f}".format( HurstEXP( mr  ) ) )
                 print ( "HurstEXP(          Trending Series ):   {0: > 12.8f}".format( HurstEXP( tr  ) ) )

                 return ( "SYNTH series demo ( on HurstEXP( ts == [ None, ] ) ) # actual numbers vary, as per np.random.randn() generator" )
            """                                                         # FIX:
            ===================================================================================================================
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :1000,QuantFX.idxH].tolist() )
            0.47537688039105963
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :101,QuantFX.idxH].tolist() )
            -0.31081076640420308
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :100,QuantFX.idxH].tolist() )
            nan
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :99,QuantFX.idxH].tolist() )

            Intel MKL ERROR: Parameter 6 was incorrect on entry to DGELSD.
            C:\Python27.anaconda\lib\site-packages\numpy\lib\polynomial.py:594: RankWarning: Polyfit may be poorly conditioned
            warnings.warn(msg, RankWarning)
            0.026867491053098096
            """
            pass;     too_short_list = 101 - len( ts )                  # MUST HAVE 101+ ELEMENTS
            if ( 0 <  too_short_list ):                                 # IF NOT:
                 ts = too_short_list * ts[:1] + ts                      #    PRE-PEND SUFFICIENT NUMBER of [ts[0],]-as-list REPLICAS TO THE LIST-HEAD
            #---------------------------------------------------------------------------------------------------------------------------
            lags = range( 2, 100 )                                                              # Create the range of lag values
            tau  = [ np.sqrt( np.std( np.subtract( ts[lag:], ts[:-lag] ) ) ) for lag in lags ]  # Calculate the array of the variances of the lagged differences
            #oly = np.polyfit( np.log( lags ), np.log( tau ), 1 )                               # Use a linear fit to estimate the Hurst Exponent
            #eturn ( 2.0 * poly[0] )                                                            # Return the Hurst exponent from the polyfit output
            """ ********************************************************************************************************************************************************************* DONE:[MS]:ISSUE / FIXED ABOVE
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ : QuantFX.aMinPTR,QuantFX.idxH] )
            C:\Python27.anaconda\lib\site-packages\numpy\core\_methods.py:82: RuntimeWarning: Degrees of freedom <= 0 for slice
              warnings.warn("Degrees of freedom <= 0 for slice", RuntimeWarning)
            C:\Python27.anaconda\lib\site-packages\numpy\core\_methods.py:94: RuntimeWarning: invalid value encountered in true_divide
              arrmean, rcount, out=arrmean, casting='unsafe', subok=False)
            C:\Python27.anaconda\lib\site-packages\numpy\core\_methods.py:114: RuntimeWarning: invalid value encountered in true_divide
              ret, rcount, out=ret, casting='unsafe', subok=False)
            QuantFX.py:23034: RuntimeWarning: divide by zero encountered in log
              return ( 2.0 * np.polyfit( np.log( lags ), np.log( tau ), 1 )[0] )                  # Return the Hurst exponent from the polyfit output ( a linear fit to estimate the Hurst Exponent )

            Intel MKL ERROR: Parameter 6 was incorrect on entry to DGELSD.
            C:\Python27.anaconda\lib\site-packages\numpy\lib\polynomial.py:594: RankWarning: Polyfit may be poorly conditioned
              warnings.warn(msg, RankWarning)
            0.028471879418359915
            |
            |
            |# DATA:
            |
            |>>> QuantFX.DATA[ : QuantFX.aMinPTR,QuantFX.idxH]
            memmap([ 1763.31005859,  1765.01000977,  1765.44995117,  1764.80004883,
                     1765.83996582,  1768.91003418,  1771.04003906,  1769.43994141,
                     1771.4699707 ,  1771.61999512,  1774.76000977,  1769.55004883,
                     1773.4699707 ,  1773.32995605,  1770.08996582,  1770.20996094,
                     1768.34997559,  1768.02001953,  1767.59997559,  1767.23999023,
                     1768.41003418,  1769.06994629,  1769.56994629,  1770.7800293 ,
                     1770.56994629,  1769.7800293 ,  1769.90002441,  1770.44995117,
                     1770.9699707 ,  1771.04003906,  1771.16003418,  1769.81005859,
                     1768.76000977,  1769.39001465,  1773.23999023,  1771.91003418,
                     1766.92004395,  1765.56994629,  1762.65002441,  1760.18005371,
                     1755.        ,  1756.67004395,  1753.48999023,  1753.7199707 ,
                     1751.92004395,  1745.44995117,  1745.44995117,  1744.54003906,
                     1744.54003906,  1744.84997559,  1744.84997559,  1744.34997559,
                     1744.34997559,  1743.75      ,  1743.75      ,  1745.23999023,
                     1745.23999023,  1745.15002441,  1745.31005859,  1745.47998047,
                     1745.47998047,  1749.06994629,  1749.06994629,  1748.29003906,
                     1748.29003906,  1747.42004395,  1747.42004395,  1746.98999023,
                     1747.61999512,  1748.79003906,  1748.79003906,  1748.38000488,
                     1748.38000488,  1744.81005859,  1744.81005859,  1736.80004883,
                     1736.80004883,  1735.43005371,  1735.43005371,  1737.9699707
                     ], dtype=float32
                    )
            |
            |
            | # CONVERTED .tolist() to avoid .memmap-type artifacts:
            |
            |>>> QuantFX.DATA[ : QuantFX.aMinPTR,QuantFX.idxH].tolist()
            [1763.31005859375, 1765.010009765625, 1765.449951171875, 1764.800048828125, 1765.8399658203125, 1768.9100341796875, 1771.0400390625, 1769.43994140625, 1771.469970703125, 1771.6199951171875, 1774.760
            859375, 1743.75, 1743.75, 1745.239990234375, 1745.239990234375, 1745.1500244140625, 1745.31005859375, 1745.47998046875, 1745.47998046875, 1749.0699462890625, 1749.0699462890625, 1748.2900390625, 174
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ : QuantFX.aMinPTR,QuantFX.idxH].tolist() )
            C:\Python27.anaconda\lib\site-packages\numpy\core\_methods.py:116: RuntimeWarning: invalid value encountered in double_scalars
              ret = ret.dtype.type(ret / rcount)

            Intel MKL ERROR: Parameter 6 was incorrect on entry to DGELSD.
            C:\Python27.anaconda\lib\site-packages\numpy\lib\polynomial.py:594: RankWarning: Polyfit may be poorly conditioned
              warnings.warn(msg, RankWarning)
            0.028471876494884543
            ===================================================================================================================
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :1000,QuantFX.idxH].tolist() )
            0.47537688039105963
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :101,QuantFX.idxH].tolist() )
            -0.31081076640420308
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :100,QuantFX.idxH].tolist() )
            nan
            |
            |>>> QuantFX.HurstEXP( QuantFX.DATA[ :99,QuantFX.idxH].tolist() )

            Intel MKL ERROR: Parameter 6 was incorrect on entry to DGELSD.
            C:\Python27.anaconda\lib\site-packages\numpy\lib\polynomial.py:594: RankWarning: Polyfit may be poorly conditioned
            warnings.warn(msg, RankWarning)
            0.026867491053098096
            """
            return ( 2.0 * np.polyfit( np.log( lags ), np.log( tau ), 1 )[0] )

有关属性的详细帖子和讨论是 here