2. 随机变量的数学统计特征

digraph {
    序列数据->自相关序列
    序列数据->序列不相关序列（白噪音）
}

import statsmodels.tsa.api as smt
"""
tsa是Time Series analysis缩写
tsa的stattools（统计工具）提供了计算acf和pacf以及后面要用到的adfuller单位根检验函数
使用help(smt.stattools.acf)可以查看相关参数设置
"""

# 计算自相关系数，这里设置滞后项为5期,默认是40期滞后
acf=smt.stattools.acf(df,nlags=5)
# 计算偏自相关系数
pacf=smt.stattools.pacf(df,nlags=5)
print(f'自相关系数为：{acf};\n偏自相关系数为：{pacf}')

>>>
自相关系数为：[1. 0.99098764 0.98189466 0.97312885 0.96252012 0.95335064]
偏自相关系数为：[ 1. 0.99231072 -0.01047826  0.01620047 -0.12635305  0.09200772]


from statsmodels.graphics import tsaplots
tsapltots.plot_acf(Train_Y,use_vlines=True,lags=80)
print("60 以后tsapltots判定为acf==0 (alpha=0.05,95%置信度)")
plt.show()

import statsmodels.tsa.api as smt

pacf=smt.stattools.pacf(df,nlags=5)
# pacf .shape ==df.shape 

from statsmodels.graphics import tsaplots
tsaplots.plot_pacf(
    x,
    ax=None,
    lags=None,
    alpha=0.05,
    method='ywunbiased',
    use_vlines=True,
    title='Partial Autocorrelation',
    zero=True,
    vlines_kwargs=None,
    **kwargs,
)

print("60 以后tsapltots判定为acf==0 (alpha=0.05,95%置信度)")
plt.show()

import statsmodels
q_statarray,p_value =statsmodels.tsa.stattools.q_stat(x, nobs, type='ljungbox')

q_statarray  # shape =x.shape
>>>

p_value  # shape =x.shape
>>>

# p_valeue 一般查看最后一个值
# 当 p_valeue > 5%(0.05) ,则接收假设，认为序列是白噪音 
# 当 p_valeue <5% , 则拒绝假设，认为序列是非白噪音 

digraph {
    序列数据->非平稳序列
    序列数据->平稳序列
    平稳序列->强平稳序列（理论上）
    平稳序列->弱平稳序列[color=red,label="一般指的是"]
}

digraph a{
    假设存在单位根->平稳序列[label="不存在 p_value<0.05",color=blue]
    假设存在单位根->非平稳序列[label="存在 p_value>0.05",color=red]
}

from arch.unitroot import ADF

ADF_instance=ADF(x,lags=10)

print (ADF_instance)
>>>
   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -8.635
P-value                         0.000
Lags                               10
-------------------------------------

Trend: Constant
Critical Values: -3.43 (1%), -2.86 (5%), -2.57 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

# 若按照5% 致信区间判定

if ADF_instance.stat<ADF_instance.critical_values["5%"]:
    print("拒绝假设，不存在单元根，序列是宽平稳的，是平稳序列")

if ADF_instance.p_value<0.05:
    print("拒绝假设，不存在单元根，序列是宽平稳的，是平稳序列")

stattools.ccovf(x, y[, adjusted, demean])

Calculate the crosscovariance between two series.

stattools.ccf(x, y[, adjusted])

The cross-correlation function.

#cointegration 协整性
coint_t,p_value,crit_value=statsmodels.tsa.stattools.coint(y0,y1,trend ='c',method ='aeg',maxlag = None,autolag ='aic',return_results = None)

if p_value<0.01:
    print("y0 y1 存在协整关系")
y0:

statsmodels.tsa.vector_ar.vecm.coint_johansen(endog, det_order, k_ar_diff)

def coint_johansen(endog, det_order, k_ar_diff):
    """
    Johansen 检验

    参数
    ----------
    endog : array_like (nobs_tot x neqs)
        Data to test
    det_order : int
        * -1 - no deterministic terms
        * 0 - constant term
        * 1 - linear trend
    k_ar_diff : int, nonnegative
        Number of lagged differences in the model.

    Returns
    -------
    result : JohansenTestResult
        An object containing the test's results. The most important attributes
        of the result class are:

        * trace_stat and trace_stat_crit_vals
        * max_eig_stat and max_eig_stat_crit_vals

    """
    import warnings
    if det_order not in [-1, 0, 1]:
        warnings.warn("Critical values are only available for a det_order of "
                      "-1, 0, or 1.", category=HypothesisTestWarning)
    if endog.shape[1] > 12:  # todo: test with a time series of 13 variables
        warnings.warn("Critical values are only available for time series "
                      "with 12 variables at most.",
                      category=HypothesisTestWarning)

    from statsmodels.regression.linear_model import OLS

    def detrend(y, order):
        if order == -1:
            return y
        return OLS(y, np.vander(np.linspace(-1, 1, len(y)),
                                order+1)).fit().resid

    def resid(y, x):
        if x.size == 0:
            return y
        r = y - np.dot(x, np.dot(np.linalg.pinv(x), y))
        return r

    endog = np.asarray(endog)
    nobs, neqs = endog.shape

    # why this?  f is detrend transformed series, det_order is detrend data
    if det_order > -1:
        f = 0
    else:
        f = det_order

    endog = detrend(endog, det_order)
    dx = np.diff(endog, 1, axis=0)
    z = lagmat(dx, k_ar_diff)
    z = z[k_ar_diff:]
    z = detrend(z, f)

    dx = dx[k_ar_diff:]

    dx = detrend(dx, f)
    r0t = resid(dx, z)
    # GH 5731, [:-0] does not work, need [:t-0]
    lx = endog[:(endog.shape[0]-k_ar_diff)]
    lx = lx[1:]
    dx = detrend(lx, f)
    rkt = resid(dx, z)  # level on lagged diffs
    # Level covariance after filtering k_ar_diff
    skk = np.dot(rkt.T, rkt) / rkt.shape[0]
    # Covariacne between filtered and unfiltered
    sk0 = np.dot(rkt.T, r0t) / rkt.shape[0]
    s00 = np.dot(r0t.T, r0t) / r0t.shape[0]
    sig = np.dot(sk0, np.dot(inv(s00), sk0.T))
    tmp = inv(skk)
    au, du = np.linalg.eig(np.dot(tmp, sig))  # au is eval, du is evec

    temp = inv(np.linalg.cholesky(np.dot(du.T, np.dot(skk, du))))
    dt = np.dot(du, temp)

    # JP: the next part can be done much  easier
    auind = np.argsort(au)
    aind = np.flipud(auind)
    a = au[aind]
    d = dt[:, aind]
    # Normalize by first non-zero element of d, usually [0, 0]
    # GH 5517
    non_zero_d = d.flat != 0
    if np.any(non_zero_d):
        d *= np.sign(d.flat[non_zero_d][0])

    #  Compute the trace and max eigenvalue statistics
    lr1 = np.zeros(neqs)
    lr2 = np.zeros(neqs)
    cvm = np.zeros((neqs, 3))
    cvt = np.zeros((neqs, 3))
    iota = np.ones(neqs)
    t, junk = rkt.shape
    for i in range(0, neqs):
        tmp = np.log(iota - a)[i:]
        lr1[i] = -t * np.sum(tmp, 0)
        lr2[i] = -t * np.log(1-a[i])
        cvm[i, :] = c_sja(neqs - i, det_order)
        cvt[i, :] = c_sjt(neqs - i, det_order)
        aind[i] = i

    return JohansenTestResult(rkt, r0t, a, d, lr1, lr2, cvt, cvm, aind)