多系列でのACF, PACFの加重平均

時系列データがはたしてどんなデータかを眺めるときの基本統計量として、自己相関(ACF, autocorrelation function), 偏自己相関(PACF, partial autocorrelation function)という指標があります。ACF、PACFについては以下がわかりやすいです。

https://sites.google.com/site/leihcrev/r/estimate-model-of-timeseries
http://st-hakky.hatenablog.com/entry/2017/05/25/123756
https://logics-of-blue.com/python-time-series-analysis/
(Logics-of-blueは「時空間分析と状態空間モデルの基礎」の著者によるblogです。本は注文中なのですがまだ入手できていません…。）
https://momonoki2017.blogspot.com/2018/03/python7.html
（対話形式でとっつきやすい）

多系列でのACF加重平均

普通時系列解析というと1か2つの時系列データを深く解析していくことが多いと思いますが、ちょっと多系列のデータを触る場面があり、ACFの加重平均をとってみました。独立な200個の仮想世界における株価の推移、みたいなイメージのデータでした。
なお、どんなデータかを見るために探索的に行なっただけで**統計的に正しいわけではないと思います。**詳しい方、コメントをいただけると嬉しいです。

acfの計算自体はstatsmodels.api.tsa.stattools.acfを使いました。この関数は基本的には一つの系列に対して用いられることが想定されているので、pd.groupbyと組み合わせます。さらに、このときのdfの各列には異なる種類のデータが入っていたので、それぞれに対してACF加重平均を求めました。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

df=pd.DataFrame(np.vstack((
    [['a', i, np.random.rand(), np.random.rand()] for i in range(100)],
    [['b', i, np.random.rand(), np.random.rand()] for i in range(100)],
    [['c', i, np.random.rand(), np.random.rand()] for i in range(100)],
    [['d', i, np.random.rand(), np.random.rand()] for i in range(100)],
    ))
    , columns=['group_id', 'time','A','B'])

df=df.set_index(['group_id', 'time']).astype(float).sort_index()
# ここまではdfの準備

grouped=df.groupby('group_id')

acf_sum_dict={}
pacf_sum_dict={}
n_k_sum_dict={}

MAX_LAG=20

for pat_id, group in grouped:
    for col_name, series in group.iteritems():
        group_acf = sm.tsa.stattools.acf(series.values, nlags=MAX_LAG)
        group_acf = group_acf[~np.isnan(group_acf)] #NaNを取り除いてから
        group_acf_filled = np.hstack([group_acf, np.zeros(MAX_LAG+1-len(group_acf))]) #0埋め

        group_pacf = sm.tsa.stattools.pacf(series.values, nlags=MAX_LAG)
        group_pacf = group_pacf[~np.isnan(group_pacf)]
        group_pacf_filled = np.hstack([group_pacf, np.zeros(MAX_LAG+1-len(group_pacf))])
                
        group_n_k = [max(0,len(series)-k) for k in range(MAX_LAG+1)] #0埋めされるように一工夫

        acf_sum_dict[col_name] = acf_sum_dict.get(col_name, np.zeros(MAX_LAG+1)) + group_acf_filled * group_n_k
        pacf_sum_dict[col_name] = pacf_sum_dict.get(col_name, np.zeros(MAX_LAG+1)) + group_pacf_filled * group_n_k
        n_k_sum_dict[col_name] = n_k_sum_dict.get(col_name, np.zeros(MAX_LAG+1)) + group_n_k

for col_name in df.columns.values:
    left = np.arange(0,MAX_LAG+1)
    acf = acf_sum_dict[col_name] / n_k_sum_dict[col_name]
    pacf = pacf_sum_dict[col_name] / n_k_sum_dict[col_name]

    plt.bar(left, acf, color='blue')
    plt.title('acf_'+str(col_name))
    plt.show()

    plt.bar(left, pacf, color='red')
    plt.title('pacf_'+str(col_name))
    plt.show()

するとこんな感じで、dfの各列に対する、groupごとのacfを加重平均したものが出てきます。
上のdfと違うときの出力なので、グラフの雰囲気はちょっと違います。各groupはランダムな長さにしたのと少ないデータ数で実験したので、グラフ右側でのデータ数が小さくなってpacfがかなりばらついています。

かなり

ニッチな場面だと思いますが、なにかのヒントになれば。