UKgasで時系列データ分析

今まで時系列のplotはしても分析らしい分析に触れる機会がなかったが、セミナーで習ったことを振り返りついでになぞる。

原理の部分は次回


ホワイトノイズ(不規則なノイズ)の生成

```python
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
import statsmodels.api as sm

a = np.random.normal(0, 1, 100)
df_whitenoise = pd.DataFrame(a, columns = ["whitenoise"])
df_whitenoise.head()
df_whitenoise.plot()
```

whitenoise
0 	-0.096781
1 	0.140867
2 	-0.624550
3 	-1.946654
4 	0.274500

![index.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/db6cfc5a-d16e-9206-d6b4-71f8371fab5c.png)

ランダムウォークの生成
ホワイトノイズの値の累積和(前回の値と足し算)することで、ランダムな時系列変化があるようなグラフを作る

```python
df_randomwalk = df_whitenoise.assign(whitenoise = np.cumsum(df_whitenoise.whitenoise))
plt.plot(df_randomwalk)
```
![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b6d86c1d-cf14-6dbe-b089-911b5b5624f9.png)

下がっていくような傾向があるように見えるが、実際は無意味に動いているだけ。
投資などの予想が難しいのはこのため。

値をズラすshift(periods=1)
値をズラしたデータを作ることで自己相関の計算がしやすくなる。

```python
df_whitenoise["shift1"] = df_whitenoise["whitenoise"].shift(periods=1)
df_whitenoise["shift2"] = df_whitenoise["whitenoise"].shift(periods=2)
df_whitenoise.head()
```

![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d6b2cf5f-f337-83ec-96bf-d7e0346d589e.png)

自己相関を見るACF
自己相関関数は、k時間単位離れた時系列の観測値間の相関（ytとyt-k）を表す

```python
sm.graphics.tsa.plot_acf(df_whitenoise.whitenoise, alpha=0.05,lags=20)
```

![06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/f1631f14-6eb4-def3-7b02-168f819d1bad.png)

最初のスパイクが優位に検知されている。
その後のスパイクは有意水準95%の範囲内に収まっている。


偏自己相関PCAF
偏自己相関関数は、遅れが短いその他のすべての項（yt-1、yt-2、...、yt-k-1）の存在を考慮して調整した後に計算される、k時間単位離れた時系列の観測値間の相関（ytとyt-k）を表す

```python
sm.graphics.tsa.plot_pacf(df_whitenoise.whitenoise, alpha=0.05,lags=20)
```

![07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/3ed80dd5-3b1b-a41a-c750-89ecdbdc6b2f.png)

移動平均過程　MA
時系列はすべての時刻で同じ分布からデータが出てくるわけではない。
時刻が異なれば、分布も異なる。
過去の値を加えていくことで現在の分布を表している。

```
df_whitenoise["ma"] = df_whitenoise["whitenoise"] + df_whitenoise["shift1"]*0.5 + df_whitenoise["shift2"]*0.8 + 1
df_whitenoise = df_whitenoise.dropna()
plt.plot(df_whitenoise[["ma"]])
```

![10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/957e5265-639c-0f7d-1a6b-acf087234134.png)

```
sm.graphics.tsa.plot_acf(df_whitenoise[["ma"]], lags=20)
sm.graphics.tsa.plot_pacf(df_whitenoise[["ma"]], lags=20)
```

![11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4c936cfa-c82b-fc08-33c6-892a37db2ab7.png)


![12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/7b08988b-c4b1-38e8-337e-b9fb8f34ac45.png)

自己回帰過程 AR
過去のデータから影響を受ける。
定数部は確定できない影響を説明している。

```python
df_whitenoise["ar"] = 0
df_whitenoise["shift1"] = 0
df_whitenoise["shift2"] = 0

for i in range(100):
    df_whitenoise.iloc[i, 1] = df_whitenoise.iloc[i, 0] + 0.8*df_whitenoise.iloc[i, 2] + 0.3*df_whitenoise.iloc[i, 3] + 2
    df_whitenoise.iloc[i, 2] = df_whitenoise.iloc[i-1, 1]
    df_whitenoise.iloc[i, 3] = df_whitenoise.iloc[i-2, 1]

plt.plot(df_whitenoise[["ar"]])
sm.graphics.tsa.plot_acf(df_whitenoise[["ar"]], lags=20)
sm.graphics.tsa.plot_pacf(df_whitenoise[["ar"]], lags=20)
```
![20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0e06f863-93ac-bd61-c975-82b99fa5ea9b.png)
ACF
![21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/c84f0368-9f74-0f14-090b-b925d023f4c5.png)
PACF
![22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/91c1c611-5675-ee51-8fbd-29d214144d59.png)


#UKgasで時系列予想
(同志社大のデータサイエンス研究室のHPでRでおなじようなことをやっている)

UKgasはRからcsvを作っておく

```R
gas_csv <- data.frame(UKgas)

year <- seq(1960,1986,1)
mon <- c(1,4,7,10)
day <- c(1)

limit <- (nrow(gas_csv)/4)
time <- NULL

for (i in 1:limit){
N1<- paste0(year[i],"/",mon[1],"/",day)
N2<- paste0(year[i],"/",mon[2],"/",day)
N3<- paste0(year[i],"/",mon[3],"/",day)
N4<- paste0(year[i],"/",mon[4],"/",day)
time <- c(time,N1,N2,N3,N4)
}

write.csv(data.frame(time=time, UKgas=gas_csv), "UKgas.csv" ,row.names=F)
```

パイソンに戻る
とりあえずプロット

```python
df = pd.read_csv("UKgas.csv", index_col=0)
df.index = pd.to_datetime(df.index)
plt.plot(df)
```
![30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/01beefbb-d42f-9314-27bb-56275cfdf728.png)

Rでやるならこんな感じ

```R
UKGAS <- data.frame(time=time, UKgas=gas_csv)
ts.plot(UKGAS[,2])
```
![31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/90c5b6e6-5a67-dab7-da5f-31b6edcf3eff.png)

ガスは四半期で周期的に動いているのがわかる(ギザギザしている)
冬は消費量が多く夏は少ない。
増え方も大きい。

対数をとって差分(足し引き)のカタチにしてわかりやすくする。

```python
df["UKgas"] = np.log(df["UKgas"])
plt.plot(df)
sm.graphics.tsa.plot_acf(df, lags=40)
sm.graphics.tsa.plot_pacf(df, lags=40)
```
![40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/13d30ba9-4f07-0148-26cd-95ed32d1388c.png)

ACF
四半期の周期性がわかる
![41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b6daf395-9b6b-96ed-c81b-c5bd368bf59c.png)

PACF
最初の方のスパイクが反応している

![42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/fa5cbfc2-4b3a-24dc-388d-e70eadcaf139.png)


```
fit = sm.tsa.SARIMAX(df, order=(5,0,0), seasonal_order=(0,0,0,4), 
                     enforce_invertibility=False,
                     enforce_stationarity=False).fit()

fit.summary()
plt.plot(fit.resid)
```

SARIMAXを調整してAICやBICが小さくなるようなモデルを見つける。
![70.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/bb7a7586-ab41-5f27-531a-ba1b15824881.png)

誤差を見ると最初は反応しているが、あとは小さい値。
1969~1977年に誤差が発生している
![60.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/761bf37a-2d32-3108-5a8f-759cd68d6979.png)

predictで予想値を二年分出す

```python
pred = fit.predict('1986-10-01', '1988-10-01')
plt.plot(df)
plt.plot(pred, "r")
```
![61.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/1a48a744-15e7-5e43-1f32-b3d5f9c2e972.png)

対数処理をnp.expで戻す

![62.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/eff375f8-0787-faed-f419-93e53ab49a45.png)


予測値が見えました。
上昇傾向にありそうです。
周期性も再現できていそうです。