LoginSignup
1

More than 3 years have passed since last update.

【要素分解入門】時系列解析の手法をRとpythonで並べてみる♬

Posted at

今回は、前回の時系列分解をさらに追及してみた。
すなわち、モデルは以下のとおり
$Y_t=T_t + S_t + e_t$
または、
$Y_t=T_tS_te_t$
として、時系列データを分解するdecompose()手法である。
前回は、それぞれのデータを同時描画していたが、今回は例えばノイズ除去した純粋な音声や振動などの再生や利用を目的に、それぞれ保存したり再生したいという動機から実施してみた。
データは参考③のサイトの死亡数の速報値と参考④の2018年までの過去データである。
今回の参考は以下のとおりである。
参考①はR言語のものであり、参考②はpythonのものである。
ここでは独立に数値化してcsvファイルに保存するまでをまとめる。
【参考】
Extracting Seasonality and Trend from Data: Decomposition Using R
TIME SERIES DECOMPOSITION & PREDICTION IN PYTHON
統計で見る日本 e-Statは、日本の統計が閲覧できる政府統計ポータルサイトです
データセット一覧(5-4死亡月別にみた年次別死亡数及び死亡率(人口千対))@統計で見る日本
なお、今回のコードとデータは以下に置きました。
MuAuan/R_Audio

やったこと

・Rでのdecompose
・pythonでのdecompose

・Rでのdecompose

前回のほぼまんまですが、decompose_IIP$seasonalという表記でseasonalのデータを指定できます。

data <- read.csv("./industrial/sibou_p.csv",encoding = "utf-8")
IIP <- ts(data,start=c(2008),frequency=12)
w_list <- c("seasonal","random","trend","observed")
decompose_IIP <- decompose(IIP)
for (i in w_list){
  if (i=="seasonal"){
    plot(decompose_IIP$seasonal)
  }else if(i=="random"){
    plot(decompose_IIP$random)
  }else if(i=="trend"){
    plot(decompose_IIP$trend)
  }else if(i=="observed"){
    plot(IIP)
  }
}
seasonal random
IIP_compose_sibou_seasonal.png IIP_compose_sibou_random.png
trend observed
IIP_compose_sibou_trend.png IIP_compose_sibou_observed.png

もう一つのコード
まず、オリジナルデータを描画します。
次に、今回のデータは12か月周期なのでトレンド曲線として12移動平均を取ります。
そして、上記のオリジナルデータに重畳してlinesを引きます。

#install.packages("forecast")
library(forecast)
plot(as.ts(IIP))

trend_beer = ma(IIP, order = 12, centre = T) #4
lines(trend_beer)

plot_line_original.png

plotはそのトレンド曲線だけを描画しています。

plot(as.ts(trend_beer))

trend.png
オリジナルデータからトレンドデータを引き算してプロットします。

detrend_beer = IIP - trend_beer
plot(as.ts(detrend_beer))

detrend.png
そして、トレンドを差し引いたノイズが乗っているデータを12個ずつのデータに区分して、それらを平均したものを季節変動として計算します。
それを13個連ねて描画しています。
【参考】
Row and Column Summaries
matrix

m_beer = t(matrix(data = detrend_beer, nrow = 12))
seasonal_beer = colMeans(m_beer, na.rm = T)
plot(as.ts(rep(seasonal_beer,13)))

seasonal_bear13.png
最後に、ランダムな変動分を求めます。
結果はdecompose()関数で求めたものと一致しました。
こんな計算で求めているんですね。

random_beer = IIP - trend_beer - seasonal_beer
plot(as.ts(random_beer))

calc_random.png
最後にいろいろ利用するために、以下のようにcsvファイルに出力します。

write.csv(as.ts(rep(seasonal_beer,1300)), file = 'decompose_seasonal1300.csv')

・pythonでのdecompose

pythonでも同じような関数が用意されています。
コードは以下のとおりです。

import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.seasonal import STL
import matplotlib.pyplot as plt

# Set figure width to 12 and height to 9
plt.rcParams['figure.figsize'] = [12, 9]
df = pd.read_csv('sibou_.csv')
series = df['Price']
print(series)

cycle, trend = sm.tsa.filters.hpfilter(series, 144)
fig, ax = plt.subplots(3,1)
ax[0].plot(series)
ax[0].set_title('Price')
ax[1].plot(trend)
ax[1].set_title('Trend')
ax[2].plot(cycle)
ax[2].set_title('Cycle')
plt.show()

Figure_1_py.png

pythonでももう一つの方法
個人的には、分解の結果は上の方法だとcycleに対してもう一段の処理が必要で、こちらはランダム分の処理までできているのでこちらの方がリーズナブルに見えます。
また、Rと比較しても、少し異なりseasonalは年々振動幅が大きくなっている様子を見せています(Rはこの変動を平均化処理したので消えているのは当たり前ですが。。。)。場合により、この手法が一番汎用性が高いと言えると思います。
【参考】
Seasonal-Trend decomposition using LOESS (STL)

stl=STL(series, period=12, robust=True)
result = stl.fit()
chart = result.plot()
chart1= result.plot(observed=False, resid=False)
chart2= result.plot(trend=False, resid=False)
plt.show()

Observed & Trend & Seasonal & Residual
Figure_2_py.png
Trend & Seasonal
Figure_3_py.png
Obeserved & Seasonal
Figure_4_py.png
また、以下のようにすると、ばらばらのグラフとして一つのグラフに出力できます。
※pl1...などの数値データを利用するため変数に代入しています

pl1 = result.observed
pl2 = result.trend
pl3 = result.seasonal
pl4 = result.resid

plt.plot(pl1)
plt.plot(pl2)
plt.plot(pl3)
plt.plot(pl4)
plt.show()
plt.close()

pl1-4.png

また、以下のようにすると簡単に出力できます。

pl5 = result.seasonal.plot()
plt.show()
plt.close()
...
seasonal resid
seasonal.png resid.png
trend observed
trend.png observed.png

そして、csvファイル出力は以下のようにやります。

import csv
with open('sample_pl2.csv', 'w') as f:
    writer = csv.writer(f, delimiter='\n')
    writer.writerow(pl2)

まとめ

・Rとpythonでdecompose()のやり方を並べてみた
・それぞれのコードは結果は少し異なるがほぼ同等のことが出来て、ノイズやトレンドが乗った周期的な変動を純粋な周期的な変動として取り出すことが出来た
・各要素分解された成分をcsvデータに書き出すことが出来、二次利用出来る
・実は今回は死亡数の経年変化を例として3月までの速報値を使ったが、その範囲では2011年のような異常値は見えていない(むしろ減少傾向)

・次に、周期的な変動要素を音に変換し聞いてみることとする
・非周期な変動要素についても(検定はともかく)やってみる

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1