#はじめに
前回に引き続き、Rのデータ解析の勉強会に出席しました
今回は時系列データです、Pythonだとどうするの？を復習を兼ねて記録に残したいと思います
かなり不満足な結果となりましたが、統計の専門ではない人間がサクッと調べた結果ということで
RとPython（の統計関係のライブラリー充実度かな？）の違いを実感しました
Rは統計関係充実していますね~
#環境
windows10 Pro 64bit
R :3.6.1
RStudio :1.2.1335
Google Colaboratoryを使用
Python :3.6.8
#何するの？
Rに実装されている？データセットUKgasを用いて時系列データ解析を行う
①差分データをプロットしてみる
②対数にしてみる
③自己相関関数（acf:autocorrelation function)で周期性を見てみる
④スペクトル分析をしてみる
⑤単位根検定をしてみる
⑥ARモデル（自己回帰モデル）を見てみる
独立性検定 Box-Pierce検定をしてみる
正規性検定 Ljung-Box検定をしてみる
予測をしてみる
⑦差分を対数にしてARモデルを試す
⑧ARMAモデル（自己回帰移動平均モデル）を試してみる
#Rのコード
データセットの読み込みと確認
#データの読み込み
data(UKgas)
class(UKgas)#データ形式の確認 ts 時系列 (time series)
[1] "ts"
UKgas#1960年から1986年までの四半期ごとのガスの使用量？
Qtr1 Qtr2 Qtr3 Qtr4
1960 160.1 129.7 84.8 120.1
1961 160.1 124.9 84.8 116.9
...
1985 1087.0 534.7 281.8 787.6
1986 1163.9 613.1 347.4 782.8
生データプロット
##①差分データをプロット
UKgas %>% diff() %>% plot()
##②対数でプロット
UKgas %>% log() %>% diff() %>% plot()
##③自己相関関数で周期性を見てみる
acf(UKgas)
1年単位で自己相関が高く、周期性がみられる
Lag0は自分自身との相関なので相関係数1
##④スペクトル分析をしてみる
1.0に周期がありますね
UKgas %>% spec.pgram()
##⑤単位根検定をしてみる
時系列が単位根を持つかどうか（ランダムウォークであるか） を調べること
単位根検定とは、「|a| = 1」であるという帰無仮説検定である。
Phillips-Perron 検定
p値が0.01なので、単位根があるという帰無仮説が棄却される
PP.test(UKgas)
Phillips-Perron Unit Root Test
data: UKgas
Dickey-Fuller = -10.126, Truncation lag parameter = 4, p-value = 0.01
Augmented Dickey-Fuller 検定
（対立仮説は「|a| < 1 (stationary)」か「|a| > 1 (explosive)」を選択できる）
stationary：p値が0.73なので棄却できない
explosive：p値が0.26なので棄却できない
あれ？困ったぞ
adf.test(UKgas)
Augmented Dickey-Fuller Test
data: UKgas
Dickey-Fuller = -1.6079, Lag order = 4, p-value = 0.7393
alternative hypothesis: stationary
adf.test(UKgas,alternative = "explosive")
Augmented Dickey-Fuller Test
data: UKgas
Dickey-Fuller = -1.6079, Lag order = 4, p-value = 0.2607
alternative hypothesis: explosive
差分データで検定してみると棄却できた！
adf.test(diff(UKgas))
data: diff(UKgas)
Dickey-Fuller = -7.1679, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary
Warning message:
In adf.test(diff(UKgas)) : p-value smaller than printed p-value
##⑥ARモデル（自己回帰モデル）を見てみる
ARモデルの係数を算出
独立性検定 Box-Pierce検定、Ljung-Box検定をしてみる
p値が0.1378なので独立性があるとは言い切れない、相関関係がありそう
正規性検定 ジャック・ベラ検定検定をしてみる
P値が7.122e-05 勉強会では正規性はありそうとのことでしたが
帰無仮説 (H0) は標本分布が正規分布に従うことなので、棄却されるんじゃないかな？確認してみます
独立性はなさそう、正規性はありそう（ないんじゃないかな？）
予測をしてみる
95%区間
結果はいまいち
ar(UKgas)
Call:
ar(x = UKgas)
Coefficients:
1 2 3 4 5 6
0.5284 -0.2951 0.5844 0.3489 0.0436 -0.2709
Order selected 6 sigma^2 estimated as 8496
#独立性の検定 Box-Pierce検定 Ljung-Box検定
Box.test(ar.UKgas$res, type = "Ljung")
Box-Ljung test
data: ar.UKgas$res
X-squared = 2.2024, df = 1, p-value = 0.1378
# 正規性検定 ジャック・ベラ検定
ar.UKgas$res[-(1:6)]%>%
jarque.bera.test()
Jarque Bera Test
data: .
X-squared = 19.1, df = 2, p-value = 7.122e-05
#予測 12 四半期 3年先まで予測
UKgas.pred <- predict(ar.UKgas, n.ahead = 12)
$pred
Qtr1 Qtr2 Qtr3 Qtr4
1987 1053.9970 600.7282 316.9538 748.8121
1988 981.4623 546.8309 308.6325 708.0975
1989 912.3671 504.7132 297.3734 674.2334
$se
Qtr1 Qtr2 Qtr3 Qtr4
1987 92.17482 104.24961 104.25994 111.21667
1988 137.88835 147.29676 147.31935 151.18307
1989 170.10928 175.66459 175.66783 178.40438
ts.plot(UKgas,UKgas.pred$pred,
gpars=list(lt=c(1,2),col=c(1,2)))
# 95%区間（±2σ）でプロット
SE1<-UKgas.pred$pred+2*UKgas.pred$se
SE2<-UKgas.pred$pred-2*UKgas.pred$se
ts.plot(UKgas,UKgas.pred$pred,
SE1,SE2,gpars=list(lt=c(1,2,3,3),col=c(1,2,4,4)))
##差分を対数にしてでARモデルを試す
独立性がありそう、相関関係がなさそう
正規性は変わらずありそう（なさそうじゃないかな？）
プロット省略
UKgas_log_diff <- UKgas %>% log() %>% diff()
ar.UKgas_log_diff <- ar(UKgas_log_diff)
#独立性検定
Box.test(ar.UKgas_log_diff$res, type = "Ljung")
Box-Ljung test
data: ar.UKgas_log_diff$res
X-squared = 7.0796, df = 1, p-value = 0.007797
#正規性検定
ar.UKgas_log_diff$res[-(1:3)]%>%
jarque.bera.test()
Jarque Bera Test
data: .
X-squared = 163.92, df = 2, p-value < 2.2e-16
#予測
UKgas_log_diff.pred
$pred
Qtr1 Qtr2 Qtr3 Qtr4
1987 0.4464858 -0.6393982 -0.5191298 0.7510022
1988 0.4517140 -0.6265824 -0.4823045 0.7013768
1989 0.4485959 -0.6073688 -0.4532143 0.6595974
$se
Qtr1 Qtr2 Qtr3 Qtr4
1987 0.1548956 0.2068423 0.2089726 0.2105422
1988 0.2449782 0.2685604 0.2726472 0.2756092
1989 0.2951754 0.3091888 0.3143251 0.3179931
# diffをcumsumで元に戻して、expすると予測値が得られる
#Pythonで頑張ってみる
RからUKgasをcsv出力
pandasでcsvを読み込むが、index (date、Qt1,Qt2,,)が抜けていたので
indexを追加
（Rから出力する際に、write.tableのrow.names=Tにしてみたけどダメだった、、？）
プロット 読み込みOK
#RでUKgasをcsv出力
write.table(UKgas, "UKgas2.csv", sep=",", row.names=T)
#Pythonのpandasで読み込み
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
UKgas = pd.read_csv('UKgas.csv')
UKgas['date']=pd.date_range('1960','1987',freq='Q')
UKgas.set_index('date', inplace=True)
#プロットしてみる
plt.plot(UKgas)
#①差分データをプロットしてみる
差分データを作って引くか .shift()
差分データをそのまま作るか .diff() どちらも可能
4半期データなので差分1～4でプロット
plt.figure(figsize=(20,5))
plt.subplot(1,4,1)
plt.plot(UKgas - UKgas.shift(1))
#plt.plot(UKgas.diff())
plt.ylim(-1000,1000)
plt.subplot(1,4,2)
plt.plot(UKgas - UKgas.shift(2))
plt.ylim(-1000,1000)
plt.subplot(1,4,3)
plt.plot(UKgas - UKgas.shift(3))
plt.ylim(-1000,1000)
plt.subplot(1,4,4)
plt.plot(UKgas - UKgas.shift(4))
plt.ylim(-1000,1000)
#②対数にしてみる
plt.figure(figsize=(20,5))
plt.subplot(1,4,1)
plt.plot(np.log(UKgas).diff(1))
plt.ylim(-2,2)
plt.subplot(1,4,2)
plt.plot(np.log(UKgas).diff(2))
plt.ylim(-2,2)
plt.subplot(1,4,3)
plt.plot(np.log(UKgas).diff(3))
plt.ylim(-2,2)
plt.subplot(1,4,4)
plt.plot(np.log(UKgas).diff(4))
plt.ylim(-2,2)
#③自己相関関数（acf:autocorrelation function)で周期性を見てみる
重回帰と同じライブラリーstatsmodelsで自己相関関数を利用することができる
プロットも同じライブラリーで描画可能
matplotlibでplotしてしまうと折れ線グラフになりみにくい、、Rは適時切り替えてくれて便利、、
barとも違うし、このグラフはなんていうグラフなんでしょうか？
結果は同じように、4半期周期で相関がみられますね
偏自己相関をとると1年後の相関が強いことがわかります
グラフの網かけは95%区間です、これは便利
import statsmodels.api as sm
#自己相関を求める
# nlag のデフォルトは40 ラグ0～40までの41個の自己相関係数を出力 0は自分自身なので常に1
UKgas_acf = sm.tsa.stattools.acf(UKgas ,nlags = 40)
UKgas_acf
array([ 1. , 0.55962589, 0.16543426, 0.5579524 , 0.90152719,
0.49563309, 0.13309608, 0.49676976, 0.80555787, 0.434424 ,
0.10824271, 0.44067708, 0.71960196, 0.37821433, 0.08281566,
0.37936599, 0.63341685, 0.31724604, 0.0517523 , 0.32556113,
0.53930818, 0.24440861, 0.0192022 , 0.26542743, 0.44460168,
0.17427765, -0.01797088, 0.20634612, 0.34597344, 0.11226117,
-0.05301067, 0.13370736, 0.24803597, 0.04672887, -0.09449297,
0.06384167, 0.16578529, -0.01332281, -0.13163598, 0.00430394,
0.08985417])
# 自己相関のグラフ 自己相関関数（ACF: autocorrelation function）
# StatsModelsのgraphics.tsa.plot_acf() で計算とプロットを同時にできる
# 偏自己相関（PACF:Pertial Autocorrelation Function）
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(UKgas, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(UKgas, lags=40, ax=ax2)
stastmodelsのseasonal_decomposeを使うと、トレンド、季節性、残差を確認できる
# statsmodel のseasonal_decomposeで俯瞰データ出力 オリジナル、トレンド、季節性、残差
res = sm.tsa.seasonal_decompose(UKgas)
# res = sm.tsa.seasonal_decompose(passengers) # 解析結果は取得済み
original = UKgas # オリジナルデータ
trend = res.trend # トレンドデータ
seasonal = res.seasonal # 季節性データ
residual = res.resid # 残差データ
plt.figure(figsize=(8, 8)) # グラフ描画枠作成、サイズ指定
# オリジナルデータのプロット
plt.subplot(411) # グラフ4行1列の1番目の位置（一番上）
plt.plot(UKgas)
plt.ylabel('Original')
# trend データのプロット
plt.subplot(412) # グラフ4行1列の2番目の位置
plt.plot(trend)
plt.ylabel('Trend')
# seasonalデータ のプロット
plt.subplot(413) # グラフ4行1列の3番目の位置
plt.plot(seasonal)
plt.ylabel('Seasonality')
# residual データのプロット
plt.subplot(414) # グラフ4行1列の4番目の位置（一番下）
plt.plot(residual)
plt.ylabel('Residuals')
plt.tight_layout() # グラフの間隔を自動調整
#④スペクトル分析をしてみる
statsモデルにピリオドグラムがあるので試してみる、、
あれ？予想と違う、、
パワースペクトル自体はFFTとかで出せそうですが、、periodogramの使い方確認します
from statsmodels.tsa.stattools import periodogram
UKgas_sp = periodogram(UKgas)
plt.plot(UKgas_sp)
#⑤単位根検定をしてみる
wikiによると単位根検定は以下のようなものがある
太字がRで検証してみた検定
・拡張ディッキー–フラー検定
・ディッキー–フラー検定
・フィリップス–ペロン検定
・KPSS検定（英語版）
・Zivot–Andrews検定（英語版）
statsmodelsにadf検定は搭載されているのですが、、エラーが、、調査中です
検索してもRの情報が多いですね、、この辺がRの強いところかな
#⑥ARモデル（自己回帰モデル）を見てみる
statsmodelsにARが入っています、そのままですがPythonの統計関係はstatsmodelsにほぼ入っていますね～
公式サイトの確認必要ですね、、
https://www.statsmodels.org/stable/index.html
モデルを学習して予測してみます
maxlagを2（周期以下）にするとうまく予測できていない、、初期から後半のバイアスがかかっているようですね
maxlagを4にするとそれなりに、、Rでの予測はあまりよくありませんでしたが、、
4以上では、パッと見た感じでは同程度でした
from statsmodels.tsa.ar_model import AR
%matplotlib inline
ar = AR(UKgas).fit(maxlag=4,ic='aic')
#予測
ar_predict = ar.predict(end=120)
#プロット
plt.plot(ar_predict,c='red',label="predict")
plt.plot(UKgas,c='black',label='orijinal')
plt.legend()
plt.show()
独立性検定 Box-Pierce検定をしてみる
正規性検定 Ljung-Box検定をしてみる
検定関係のいい関数がpythonで見つけられませんでした
実装しているひとは、ちらほらいるんですが、少し主旨とはずれるので、、
見つけたら追記します
事前調査はRで、実装はpythonってパターンが多いんですかね？
#⑧差分を対数にしてARモデルを試す
検定した結果diffをが良好としているのですが、検定できていないので、、
diffデータを使うだけれであれば⑦と同じです
##⑨ARMAモデル（自己回帰移動平均モデル）を試してみる
MA 移動平均モデルとARの二つを足したモデル
使用するにあたっては、ARとMRの次数を決める必要があります
aicが最小になる、AR 3,MA 2の組み合わせがいい感じ
この先でうまく予想ができない、、予測値がすべてNaNに、、調べて追記します
from statsmodels.tsa import stattools as st
# ARMAモデルの次数を決める
st.arma_order_select_ic(UKgas, ic='aic', trend = 'nc')
{'aic': 0 1 2
0 NaN 1495.582257 1435.927808
1 1479.809575 1448.623957 1424.508152
2 1481.530412 1438.384552 1388.008218
3 1208.282570 1152.034673 1118.209335
4 NaN 1122.133937 1172.926357, 'aic_min_order': (3, 2)}
#ARMA予測
from statsmodels.tsa.arima_model import ARMA
arma_UKgas = ARMA(UKgas,order=[3,2]).fit(maxlags=2, ic='aic')
arma_UKgas_predict = arma_UKgas.predict(end=120)
plt.plot(arma_UKgas_predict)
#以下も追記予定、、
統計関係はRのほうが充実していることが身に染みてわかりました、、、
⑧ARIMAモデル（自己回帰和分移動平均モデル）
⑨SARIMAモデル（季節的自己回帰和分移動平均モデル