時系列解析 -自己回帰型モデル・状態空間モデル・異常検知- 1章 を読んだ
時系列解析 -自己回帰型モデル・状態空間モデル・異常検知の1章、時系列データの記述・処理を読んだので、プログラムを再現しつつまとめてみました。
2020/3/15 追記 輪読会で使用したスライドを追加
www.slideshare.net
python のデータはこちら。
季節変動の除外
まず、下記のサイトにある、アイスクリームの金額のデータ、icecream.csvをダウンロードします。 oku.edu.mie-u.ac.jp
データを読込んでみましょう。
このとき、read_csv に encoding="shift-jis" を付けないと、読込エラーが生じました。
また、これは自分の環境だけかもしれませんが、 read_csv に engine="python" を付けないと、読み込みエラーが生じました。
df_ice = pd.read_csv('【データ保存場所のpath】icecream.csv', encoding="shift-jis",engine="python") df_ice.head(10)
表示をみると、インデックスがおかしくなってますので、index(inplace = True)でインデックスを解除してやります。
# インデックス解除 df_ice.reset_index(inplace = True)
表のカラム名もおかしくなっているので、修正して、余計な文字が入っているindex 0 を削除してやります。
# カラム名をチェック df_ice.columns # 先ほど調べたカラム名を、年,月,アイスクリームに変更 df_ice = df_ice.rename(columns={'level_0':"年",'level_1':"月", '# 家計調査 家計収支編 二人以上の世帯 全国 アイスクリーム・シャーベット【円】':"アイスクリーム"}) # 余計なindex 0 の削除 df_ice = df_ice.drop(index=0) df_ice
実はこのデータ、アイスクリームの金額がオブジェクト'dype(O)'になっていますので、このままではグラフ表示も計算もできません。
数値データに修正してあげます。
# 型をfloatに変更 df_ice["アイスクリーム"] = df_ice["アイスクリーム"].astype(float) # グラフ表示 df_ice["アイスクリーム"].plot()
これでデータの前処理は完了しました。
下記のグラフを見るとわかるように、夏の金額が大きく、それ以外の金額が小さくなっており、季節的変動が非常に大きくなっています。
そこで、今回のメインテーマである、季節的変動をデータから除外してみましょう。
季節変動を除外するために、rolling関数を使い移動平均を算定します。
window=12 を指定することで、1年周期(12か月)の移動平均を算定することができます。
.shift(-6)を使っているのは、最初の年の1か月目～6か月目は前後12か月分のデータがなく、NaNになってしまうので、グラフのスタート時点をずらすためです。
df_ma = df_ice["アイスクリーム"].rolling(window=12).mean().shift(-6)
また、厳密には、周期が12か月と偶数の場合、1月～12月の移動平均は6.5か月時点に対応するので、7月時点に対応するデータを算定するためには、6.5月時点と7.5月時点に対応する二点の平均を取る必要があるので、その処理を行います。
これにより、中心化移動平均(cma Center Moving Avarage)を求めることができます。
df_cma = df_ma.rolling(window=2).mean() # グラフを重ねて表示 df_ice["アイスクリーム"].plot() df_cma.plot(color="red",style='--')
青いのが元のデータ、赤いのが中心化移動平均後のデータです。
データを乗法モデルで考えると、T（Trend 傾向変動） × S（Seasonal 季節変動） × I（Irregular）不規則変動）と考えられます。
季節変動を除外するためには、季節的変動を表す季節指数を計算で求め、元データを季節指数で割ってやればよいのです。
次のように計算します（ここはテキストにハッキリ書いてない部分も含んでいますので、間違ってたらごめんなさい）。
① A ：もともとのアイスクリームの金額 T * S * I
② B ：中心化移動平均後のアイスクリームの金額 T
③ C A/B ：(T * S * I)/ T = S * I
④ D :12か月ごとの季節指数を求める S
→月ごとの平均値を求めてその和が1200となるようにする
⑤ E A/D =(T * S * I)/S = T * I
①、②のデータはすでに用意していますので、③の元データのアイスクリームの金額を、中心化移動平均後で割ってやります。
df_orig_div_cma = df_ice["アイスクリーム"] / df_cma
季節指数を算定するための変数を、zeros(12)を使って12個のゼロで初期化した配列用意します。
s_index が季節指数、counterが季節性のない数字(おそらくI)に該当するようです。
s_index = np.zeros(12) s_index counter = np.zeros(12, dtype='i') counter
では季節指数を算定してみましょう。 季節指数は、100を1年の平均値と考えた場合の、各月の季節効果を含んだ相対的な値です。
for idx in range(len(orig_div_cma)//12): # 12か月ごとにデータを抽出 cut_orig_div_cma = orig_div_cma[idx*12:(idx+1)*12] mask = cut_orig_div_cma!= cut_orig_div_cma # numpy.whereを使用して非数(nan)を0にして加算する counter += np.where(mask, 0 , 1) s_index += np.where(mask, 0, cut_orig_div_cma) # 加算結果の各月平均 s_index /= counter # 全体を1200に合わせ季節指数を計算 s_index = s_index / s_index.sum() * 1200 s_index
季節指数
12か月分の季節関数を各年に入力してやり、季節指数で割って100をかけてやると、季節調整済みのアイスクリームの金額を求めることができます。
# 12か月分の季節関数を各年に入力 tiled_s_index = np.tile(s_index, len(orig_div_cma)//12) # 季節調整済みのデータの計算 df_adjusted_series = df_ice["アイスクリーム"] / tiled_s_index * 100
何の計算をしているのかわかりにくいかと思いますので、さきほどの計算をスプレッドシートで実施してみましたのでご参考ください。
これで季節調整済みのデータを算定することができました。
なぜこんなややこしい計算をするのでしょう？
時系列データに季節変動がある場合、成長率を見るためには前月/前期比伸び率ではなく、前年同月/前年同期比の伸び率を見る必要があります。
例えば、アイスクリーム2019年7月の売上と、2019年6月の売上を比較しても、季節的変動の影響があるため、比較しても仕方ない可能性が高いのです。
したがって、2019年7月の売上は、2018年7月の売上と比較する必要があります。
しかし、変化の激しい現代においては、前年同月比では効果が乏しい可能性もあります。
そこで、季節的変動を取り除いた数値をベースに、前月/前期比伸び率を計算するのです。
また、季節指数を使うと、年間の目標数値から、簡単に各月の目標数値を算定することができます。
以上になります、最後までお読みいただきありがとうございました。
