statsmodelsのpacf (method='yw') が遅かったから，数倍だけ早くしてみた

# tl;dr

- statsmodelsで計算するlagの数だけyule_walkerを呼び直してるので無駄があるよ
- 1回で最長ラグ分の自己相関関数とtoeplitz行列を計算するようにしたら，ちょびっとだけ早くなったよ
- yule walker方程式を解くところはどうしようもないので，劇的には早くならないよ
    - 同モジュール内に実装済のlevinson durbinのアルゴリズムに余裕で負けるよ…！

# 注意
- statsmodelsが無駄な処理をしてるのは多分，yule_walkerのアルゴリズムに修正が入った時にまとめて全体に処理を波及させるための工夫だと思う．
- 逆に，以下のコードはyule_walker関数の処理の途中の項を使うので，修正が入る度に手を入れないといけない．  
  つまりメンテナンス性は最悪．
- しかも，どうせ速さを求めるならmethod='ld'でpacfを呼び出した方がよっぽど早い．

# 該当箇所

yule_walkerを呼び出しているところ．

```python:statsmodels.tsa.stattools
def pacf_yw(x, nlags=40, method='unbiased'):
    '''(略)'''
    pacf = [1.]
    for k in range(1, nlags + 1):
        pacf.append(yule_walker(x, k, method=method)[0][-1])
    return np.array(pacf)
```

`yule_walker`は`statsmodels.regression.linear_model`の中の関数．

```python
def yule_walker(X, order=1, method="unbiased", df=None, inv=False,
                demean=True):
    """(略)"""

    method = str(method).lower()
    if method not in ["unbiased", "mle"]:
        raise ValueError("ACF estimation method must be 'unbiased' or 'MLE'")
    X = np.array(X, dtype=np.float64)
    if demean:
        X -= X.mean()                  # automatically demean's X
    n = df or X.shape[0]

    if method == "unbiased":        # this is df_resid ie., n - p
        denom = lambda k: n - k
    else:
        denom = lambda k: n
    if X.ndim > 1 and X.shape[1] != 1:
        raise ValueError("expecting a vector to estimate AR parameters")
    r = np.zeros(order+1, np.float64)
    r[0] = (X**2).sum() / denom(0)
    for k in range(1, order+1):
        r[k] = (X[0:-k] * X[k:]).sum() / denom(k)
    R = toeplitz(r[:-1])

    rho = np.linalg.solve(R, r[1:])
    sigmasq = r[0] - (r[1:]*rho).sum()
    if inv:
        return rho, np.sqrt(sigmasq), np.linalg.inv(R)
    else:
        return rho, np.sqrt(sigmasq)
```

# 高速化

ここで，`r`は単なる自己相関関数なので，`yule_walker`関数の第二引数(`order`)に影響を受けないし，
`R`は`r[k]`を行列の$(i,j)$成分について$|i-j|=k$なる要素に並べた[toeplitz行列](https://ja.wikipedia.org/wiki/%E3%83%86%E3%83%97%E3%83%AA%E3%83%83%E3%83%84%E8%A1%8C%E5%88%97)であるので，やはり`order`に影響を受けない．
すなわち，関数`yule_walker`内の`R`はある時系列に対して1度だけ計算すれば良い．

そのように改変したのが以下である．
上記のように計算をサボるためには，yule_walkerの計算をバラしてpacf_ywに入れざるをえない．
また，stats.acfからもコードを持ってきて自己相関の計算を行なっている．

```python
def mypacf(x, nlags=40, method='unbiased'):
    method = str(method).lower()
    if method not in ["unbiased", "mle"]:
        raise ValueError("ACF estimation method must be 'unbiased' or 'MLE'")
    X = np.array(x, dtype=np.float64)
    # order = nlagsのyule walkerを解く
    if X.ndim > 1 and X.shape[1] != 1:
        raise ValueError("expecting a vector to estimate AR parameters")
    r = stats.acovf(X, unbiased=(method == 'unbiased'), demean=True)
    R = toeplitz(r[:nlags])

    # 最大のlagであるnlagsまで1度だけ計算したRとrを，再利用してpacfを算出する
    # Rと違ってpythonではlag=0から返却するので，1.0を先頭につける
    return np.array([1.] + [np.linalg.solve(R[:k,:k], r[1:k+1])[-1] for k in range(1,nlags+1)])
```

# 評価

実は，最も時間がかかるのは連立方程式を解く部分なのだが，これは`R`のサイズを増やしながら順次解かざるを得ないため，今回の方法で高速化することができない．
そのため，計算速度は劇的には改善せず，実行に要する時間は高々数分の1にしかならない．

```python
import timeit
import statsmodels.tsa.arima_process
import io

# setups
stats_setup = 'import statsmodels.tsa.stattools as stats'
my_setup = '''
import numpy as np
from scipy.linalg import toeplitz
import statsmodels.tsa.stattools as stats

def mypacf(x, nlags=40, method='unbiased'):
    method = str(method).lower()
    if method not in ["unbiased", "mle"]:
        raise ValueError("ACF estimation method must be 'unbiased' or 'MLE'")
    X = np.array(x, dtype=np.float64)
    if X.ndim > 1 and X.shape[1] != 1:
        raise ValueError("expecting a vector to estimate AR parameters")
    r = stats.acovf(X, unbiased=(method=='unbiased'), demean=True)

    R = toeplitz(r[:nlags])

    return np.array([1.] + [np.linalg.solve(R[:k,:k], r[1:k+1])[-1] for k in range(1,nlags+1)])
'''

# statements
my_stmt    = 'mypacf(data, nlags={N})'
stats_stmt = 'stats.pacf(data, nlags={N})'
ld_stmt    = "stats.pacf(data, nlags={N}, method='ld')"  # levinson durbin

# the number of samples
sample = 10

for n in [100,1000,]:  # 10000は遅すぎて力尽きた．
    data = statsmodels.tsa.arima_process.arma_generate_sample([0.5,-0.5],[0.5,0.5],n)

    my_result = timeit.timeit(my_stmt.format(N=n//2), my_setup.format(n=n), number=sample, globals=globals())
    stats_result = timeit.timeit(stats_stmt.format(N=n//2), stats_setup.format(n=n), number=sample, globals=globals())
    ld_result = timeit.timeit(ld_stmt.format(N=n//2), stats_setup.format(n=n), number=sample, globals=globals())
    print(
        f'{n}',
        f'{my_result/sample:.4f}',
        f'{stats_result/sample:.4f}',
        f'{ld_result/sample:.4f}'
    )
```

サンプルが少ないので結果は多少ブレるけど，Let's note (CF-SX4) で走らせてみた感じだとこんな結果になった．

   n | my pacf | pacf_yw | levinson durbin
----:|--------:|--------:|--------:
 100 | 0.0025  | 0.0123  | 0.0012
1000 | 0.5264  | 1.3742  | 0.1021