IT業務効率化

plot_pacfで1以上、−1以下の値が出てきたのでその意味を考える

結論

Yule-Waker方程式を利用するとき、有限データのバイアスを考慮した計算をされている。

Pythonのstatsmodelsのplot_pacfはデフォルトでその計算方法を利用しており、1以上−1以下になるケースは存在する。決して誤った結果ではない。

どうしてこの現象を調べようと思ったのか?

https://community.qlik.com/t5/Qlik-Server-Side-Extensions-Documents/AirPassengers-csv/ta-p/1476807

上のリンクにあるAirPassengers.csvというデータを利用ししてPythonのstatsmodels.api..graphics.tsa.plot_pacfを使った時に、36という値が出てきてしまいました。

これは正しい値なのでしょうか?自己相関関数は1〜−1までですが、偏自己相関関数はどうだったのでしょうか?式から確認してみるか、それとも別の方法で考えるか、どちらにしろデータとして自然なものなのか疑わしいので、どうしてこのようなデータがプロットさせるのか確認したいと思います。

そもそもplot_pacfでプロットされる偏自己相関関数とは?

偏自己相関関数はデータがそのデータのラグを取った時の相関を、その間にあるラグの相関による影響を全て排除して相関を見る方法です。

https://stats.biopapyrus.jp/time-series/acf.html

こちらの現象に関してはmatplotlibのQAや、質問箱でも投げられている内容で、こちらを確認しながら進めていきたいと思います。

https://answers.splunk.com/answers/655591/why-is-the-pacf-out-of-limits-in-the-splunk-machin.html

https://jp.mathworks.com/matlabcentral/answers/122337-wrong-function-partial-autocorrelation-pacf-parcorr-greater-than-1-1

上記のやり取りからいくつかのアドバイスを抜粋しました。

  1. データ数が不十分である可能性が高い
  2. PythonとRのアルゴリズムは異なるが、どちらも間違っているわけではない
  3. Yule-Walker方程式は定常性を仮定した計算方法であるため、データから定常性が見受けられない場合は、正しい計算が行えない。

これらについて、それぞれ以下のような検証を行う。

  1. データ数を増やしてみて、絶対値1以上の値が出現しないことを確認する
  2. RとPythonで何が異なっていたのか確認する
  3. Yule-Walker方程式以外での計算方法でどのような結果が出るのか確認する。また定常性を担保するホワイトノイズではどのような絶対値1以上の結果が出ないことを確認する。

データ数を増やしてみる

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
x = np.linspace(0, 1000*np.pi, 100000*np.pi)
y1 = np.sin(x)
fig1, ax1 = plt.subplots(figsize=(13, 2.5))
fig2, ax2 = plt.subplots(figsize=(13, 2.5))
sm.graphics.tsa.plot_acf(y1, lags=50*np.pi, ax=ax1);
sm.graphics.tsa.plot_pacf(y1, lags=50*np.pi, ax=ax2, method="ywunbiased");

上記のようなコードで正弦波を作成し、それに対して偏自己相関関数(および自己相関関数)をプロットしてみた。

正弦波のデータにおいても、絶対値1以上の値が出現してしまう。このデータを増やしたらどうなるか。xの値を10倍に増やしてみた。

x = np.linspace(0, 100*np.pi, 1000*np.pi)

自己相関関数の信頼区間の幅が広がり、有意性を持つデータが増えました。しかし偏自己相関関数は側は結局絶対値1以上の値が存在し続けることに。

データ量に意味はなさそう

Rでは絶対値1以上が出ない理由

先に結論をあげると計算方法が異なるからです。また、Python側でplot_pacfの引数のmethodの値を指定することでRと同じように絶対値1以上を出現させないようにもできました。

いかにはRのコードとグラフを記載しておきます。

%load_ext rpy2.ipython
%R require(graphics)
%R -i x
%R -i y1
#%R -i AirPassengersData
%R acf(y1, plot = TRUE, lag.max=160)
%R pacf(y1, plot = TRUE, lag.max=160)

また、両方自己相関関数は同じ形をしているものの、偏自己相関関数は全く異なる形をしており、Rではひとつ前のラグに対して負の相関を持ち、それ以外は相関0となっています。一方絶対値1以上の値を持ってしまうもののPython側の計算方法では周期に合わせて相関が変化しています。(※Python側R側と書いていますが、異なっているのは計算方式です。事項でより掘り下げます。)

PythonとRで結果が異なる理由は、偏自己相関関数の計算方法が異なるから

Yule-Walker方程式以外の計算方法

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
x = np.linspace(0, 10*np.pi, 100*np.pi)
y1 = np.sin(x)
fig1, ax1 = plt.subplots(figsize=(13, 2.5))
fig2, ax2 = plt.subplots(figsize=(13, 2.5))
fig3, ax3 = plt.subplots(figsize=(13, 2.5))
sm.graphics.tsa.plot_pacf(y1, lags=50*np.pi, ax=ax1, method="ywunbiased");
sm.graphics.tsa.plot_pacf(y1, lags=50*np.pi, ax=ax2, method="ywmle");
sm.graphics.tsa.plot_pacf(y1, lags=50*np.pi, ax=ax3, method="ols");
  • yw or ywunbiased : yule walker with bias correction in denominator for acovf. Default.
  • ywm or ywmle : yule walker without bias correction
  • ols – regression of time series on lags of it and on constant
  • ld or ldunbiased : Levinson-Durbin recursion with bias correction
  • ldb or ldbiased : Levinson-Durbin recursion without bias correction

ywunbiasedがデフォルトです。これはYule-Walker方程式を利用していますが、分母にバイアス補正をしているようです。

ソースコードを確認すると全体のデータ数nにラグ数kを引いた値を分母とすると、バイアスのないYule-Walker方程式が計算できる仕組みのようです。

    for k in range(1, order+1):
        r[k] = (x[0:-k] * x[k:]).sum() / (n - k * adj_needed)
    R = toeplitz(r[:-1])
    elif method in ('yw', 'ywu', 'ywunbiased', 'yw_unbiased'):
        ret = pacf_yw(x, nlags=nlags, method='unbiased')
    elif method in ('ywm', 'ywmle', 'yw_mle'):
        ret = pacf_yw(x, nlags=nlags, method='mle')

ywmleがYule-Walker方程式そのものです。上のコードでadj_neededがFalseとなり分母がnとなります。

https://github.com/statsmodels/statsmodels/blob/65af83bb4ebb18a0e87bd4dfe959583665fcef8b/statsmodels/tsa/stattools.py#L862

https://github.com/statsmodels/statsmodels/blob/65af83bb4ebb18a0e87bd4dfe959583665fcef8b/statsmodels/regression/linear_model.py#L1322

olsはOrdinary Least Squares、最小二乗法のことのようです。全ての値は1から−1の間に納まり、共分散定常性を持つ時系列データのにみに有効であるそうです、

def pacf_ols(x, nlags=40, efficient=True, unbiased=False):
    """
    Calculate partial autocorrelations via OLS.
    Parameters
    ----------
    x : array_like
        Observations of time series for which pacf is calculated.
    nlags : int
        Number of lags for which pacf is returned.  Lag 0 is not returned.
    efficient : bool, optional
        If true, uses the maximum number of available observations to compute
        each partial autocorrelation. If not, uses the same number of
        observations to compute all pacf values.
    unbiased : bool, optional
        Adjust each partial autocorrelation by n / (n - lag).
    Returns
    -------
    ndarray
        The partial autocorrelations, (maxlag,) array corresponding to lags
        0, 1, ..., maxlag.
    See Also
    --------
    statsmodels.tsa.stattools.pacf
        Partial autocorrelation estimation.
    statsmodels.tsa.stattools.pacf_yw
         Partial autocorrelation estimation using Yule-Walker.
    statsmodels.tsa.stattools.pacf_burg
        Partial autocorrelation estimation using Burg's method.
    Notes
    -----
    This solves a separate OLS estimation for each desired lag using method in
    _. Setting efficient to True has two effects. First, it uses
    `nobs - lag` observations of estimate each pacf.  Second, it re-estimates
    the mean in each regression. If efficient is False, then the data are first
    demeaned, and then `nobs - maxlag` observations are used to estimate each
    partial autocorrelation.
    The inefficient estimator appears to have better finite sample properties.
    This option should only be used in time series that are covariance
    stationary.
    OLS estimation of the pacf does not guarantee that all pacf values are
    between -1 and 1.
    References
    ----------
    ..  Box, G. E., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015).
       Time series analysis: forecasting and control. John Wiley & Sons, p. 66
    """
    x = array_like(x, 'x')
    nlags = int_like(nlags, 'nlags')
    efficient = bool_like(efficient, 'efficient')
    unbiased = bool_like(unbiased, 'unbiased')

    pacf = np.empty(nlags + 1)
    pacf[0] = 1.0
    if efficient:
        xlags, x0 = lagmat(x, nlags, original='sep')
        xlags = add_constant(xlags)
        for k in range(1, nlags + 1):
            params = lstsq(xlags[k:, :k + 1], x0[k:], rcond=None)[0]
            pacf[k] = params[-1]
    else:
        x = x - np.mean(x)
        # Create a single set of lags for multivariate OLS
        xlags, x0 = lagmat(x, nlags, original='sep', trim='both')
        for k in range(1, nlags + 1):
            params = lstsq(xlags[:, :k], x0, rcond=None)[0]
            # Last coefficient corresponds to PACF value (see )
            pacf[k] = params[-1]

    if unbiased:
        n = len(x)
        pacf *= n / (n - np.arange(nlags + 1))

    return pacf

biasって何なのか?

https://www.sciencedirect.com/science/article/pii/S1474667016393673

有限のデータを扱う際にYule-Walker方程式に設けるバイアスのようです。

こちらの論文を読めばより理解できそうですが、それだけでブログには収まりきらない内容になりそうなので、割愛します。

Pythonのstatsmodelsのplot_pacfはデフォルトでその計算方法を利用しており、1以上−1以下になるケースは存在する。決して誤った結果ではない。

ABOUT ME
hirayuki
今年で社会人3年目になります。 日々体当たりで仕事を覚えています。 テーマはIT・教育です。 少しでも技術に親しんでもらえるよう、noteで4コマ漫画も書いています。 https://note.mu/hirayuki