【投資・R】統計や機械学習使って"金持ち"になりたい~周期分析~【第3話】

**前回までのあらすじ**

・株のデータをちょこっとあつめた
・plotした
・前処理した

**今回の目的**
・周期を知ってパターンを見つける
・パターンから今後上がるのか下がるのかを掴みたい

#周期の分析
データは前回までのスクリプトをそのまま受け継いでいます。

前回の欠損を埋めたデータをもう一度plotしてみます。

```r
plot(fill_data$ind,fill_data$move, type="l")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/3d35d1f8-6b32-53cc-0830-abaf4dcbe1ec.png)


そもそも周期性があるのかを調べてみます。
周期性はフーリエ変換という便利な数学がありまして、
どの間隔で周期があるのかを見つけてくれる関数があります。

```r
spec.pgram(fill_data$move,log="no")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/96b49ca5-ab32-e8f5-0315-dcb414d6cc7a.png)

spec.pgram関数便利。
fftという高速フーリエ変換を行ってくれる関数もありまして、
その関数からデータ数や周波数を使って計算したら同じ図を作図できます。
でもこっちのほうが楽なので。

図の確認方法ですが、データ全体に対して周期性を探したところ、ほとんど0付近に高いピークが現れています。

このピークの出る位置が周期を表してくれるのですが、今回の図からわかるのは周期が無いという事です。
もしもsinカーブのように周期が綺麗にあれば、全体の何割くらいの位置で一周期になっているのか、という事がこの図でわかるようになります。
(相変わらず説明が下手)

#もっと周期成分を深堀

今回はデータに対して直接周期を求めました。
ですが、ファイナンス系のデータは単位根過程という、過去の値±ランダム値から成り立っていることがほとんどのようです。
過去の値というのはつまりオフセットのようなものです。
本来はランダムな値がどんな周期で足しあわされていくのかを知りたいわけです。

もうすこし説明すると、
時系列とはトレンド成分と周期成分から成るものであり、
トレンド成分は過去からの積み上げで、波がどの方向に動くのかを表しています。
データ全体を線形回帰したものと考えると分かりやすいかと。
周期成分は、トレンド成分に足し合わされるパターンを持った波の事。
完全にパターンが無いノイズもある。

トレンドと周期を足したら時系列をモデル化できるので、周期やトレンドが知りたい。

これを知るために階差系列を作ります。

```r
diff_series<-fill_data$move[2:length(fill_data$move)]- fill_data$move[1:length(fill_data$move)-1]
plot(diff_series,type="l")
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/a3d5d95b-324e-66b0-f3d1-355f04744047.png)


実はもっと簡単に求めるdiff関数も用意されています。

```r
plot(diff(fill_data$move),type="l")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/8363da7a-4a58-93d5-0ffc-2976c3fcad83.png)


これでトレンド成分がなくなり、周期的に表れる成分だけになりました。
この成分に周期があるのか確認します。


このdiffからも周期成分は見つけられないのでしょうか？

```r
library(animation)
saveGIF(
for(i in seq(2,730,10)){
plot(diff(log(fill_data$move),lag=i),main=paste0("lag is ",i),type="l")
}
, movie.name="diff.gif", interval=0.8)
```


![diff.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5258599a-e18e-e446-7738-a3ed581d3e7c.gif)


ちょっと分からない。。。

ちなみにsin(x)のようなきれいな周期を持っているデータについて差分を取るとdiffのplotは元のsin(x)とそんなに変わらないまま維持される。
差分を取るとここまでグネグネと形が変化するのであれば綺麗な周期は期待できないかも。

#周期の確認

```r
spec.pgram(diff_series,log="no")
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5e9db6a0-9e64-bdc4-8cf7-11cbb2604d11.png)

周期らしきものが少しだけ見えてきました。
0.15付近に小さい山が見えます。
これはデータ全体(約3000データ)に対して0.15あたりの位置に1周期があるという事です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/45f70f7e-7252-bf2e-fb4c-2512118c7ae4.png)

0.15あたりって何データ目かと確認すると、0.12で374という値が出ています。
このデータは株式のデータであり、実はこの企業は一年のうちに上がるタイミングが定期的にありそうだ、という事は事前に考えていました。
なのでこの374はおそらく365日周期を見つけてくれたのではないかと思います。
予想していた仮説を**"それらしい数字"**で確かめることができました。

もちろん、「1年で周期があるだろう」と考えていたために374が365とむすびついてしまった感は否めません。
0.13~0.15までも念のために確認すると467日になっています。
365とは少し離れています。

ここから新しい知見が得られました。
**「実は365日周期でなく、400~460日周期なのでは??」**
どのくらい役に立つかはわかりませんが、かんがえてもみなかった新しい仮説が得られたというのは大きい収穫です。

#もっともっと深堀

階差系列のplotを見てみると、0~500あたりは値が波が大きくないのに、500~3000までにじわじわと波が大きくなっていっているように感じます。

ファイナンス系の時系列データにはよくある傾向らしいのですが、
高い値をとるようになるとばらつきが大きくなってきて、低い値では安定する。という傾向があるようです。
階差系列のplotの波が大きくなっているのは、トレンド成分が増加しているために、
"よくある傾向"が発生していると考えられます。
ばらつきの大きさの変化による影響から、周期成分がうまく抜き出せなくなります。
ばらつきの大きさを軽減するため、元のデータを対数化してばらつきを小さくします。

```r
log_move<-log(fill_data$move)
plot(log_move,type="l")
spec.pgram(log_move,log="no")
plot(diff(log_move),type="l")
spec.pgram(diff(log_move),log="no")
```



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/c0af9eb0-ffd4-e06b-03f9-f18e1ceaa9c3.png)



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/bac9f466-cba6-f7d4-c0b8-dc84b2e14b11.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/215ca5d2-d94c-0f42-16fc-38aadac8ab7b.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4653d3f9-3a0d-7745-5a47-acd528d525c3.png)



対数化しても定常性は見つけられなかったような・・・
なにか間違えているような気もしますが、このまま進めます。

#関数で簡単に分析

今までやってきたことはstl関数で簡単に実行できてしまいます。


```r
move_ts<-ts(as.numeric(fill_data$move),frequency=2)
plot(stl(move_ts,s.window="per"))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/11b89882-fcc8-4afe-d6cc-4c76ba5b4a74.png)

元データ
季節性：周期的な影響
トレンド：増加、減少の傾向
残差：季節性でもトレンドでも説明できない影響

今回frequencyを2に設定してあるため、2つのデータを一塊にしてデータの傾向を掴んでほしいと指定しています。
2つのデータから季節性を抽出できるはずもないです。
例えば四季のある1年分の気温の時系列から、季節性を抽出したいのであれば、365日を暑いと寒いの2分割に分けるため180に設定するか、
3期や4期に分けるため、120や90に分けるなどです。

2つのデータを一塊にしているため、trendの図でデータのほとんどを説明できてしまい、残差はほとんど無くなっています。
最も適切に周期を掴めた!と言えるのはどんな時かを考えてみると、
・季節性の項目で周期が確認できる
・frequencyが小さすぎない
・トレンドがなめらか
・残差が小さい
を満たしている時だと思います。

まずはfrequencyを変化させることによって、stlがどう動くのか確かめるため、
図をgifにして動かしてみます。


```r
library(animation)

saveGIF(
for(i in seq(2,730,10)){
move_ts<-ts(as.numeric(fill_data$move),frequency=i)
plot(stl(move_ts,s.window="per"),main=paste0("now " ,i))
}
, movie.name="stl.gif", interval=0.8)
```

![sample.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/ce1d1236-d9ef-99ca-ef38-9589ecf480ec.gif)

こうしてgifで確認してみるとfrequencyが300を超えたあたりからトレンド成分がなだらかになってきます。
500まで行くと直線っぽくなってきます。
行き過ぎでしょうか。
残差の成分はある程度波を打っています。




```r
i=360
move_ts<-ts(as.numeric(fill_data$move),frequency=i)
res_move_ts<-stl(move_ts,s.window="per")
plot(res_move_ts$time.series[,3])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/08042659-9697-169d-4af6-2e24b5fd9e50.png)

この残差に周期を感じてしまいます。
念のため残差成分を取り出して確認してみます。


```r
saveGIF(
for(i in seq(2,730,10)){
move_ts<-ts(as.numeric(fill_data$move),frequency=i)
res_move_ts<-stl(move_ts,s.window="per")
spec.pgram(res_move_ts$time.series[,3],log="no")
}
, movie.name="reminder.gif", interval=0.8)
```

![reminder.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/619e7352-8880-5c0c-db11-64a433620e6c.gif)



残差成分はどれだけ動かしても周期はなさそうです。







#差分に定常性などなかった

差分をとって周期をさがしたり、
stlで成分の分解を行いながら周期を見つけようとがんばりましたが、
周期は無さそうな気がしてきました。

単位根過程の検定を行ってみると

```r
library(tseries)
adf.test(fill_data$move)
```


```r
> adf.test(fill_data$move) 

        Augmented Dickey-Fuller Test

data:  fill_data$move
Dickey-Fuller = -3.1431, Lag order = 14, p-value = 0.0977
alternative hypothesis: stationary
```

p値がほとんど0です。
差分を見てみます。

```r
adf.test(diff(fill_data$move))
```

```r
> adf.test(diff(fill_data$move))

        Augmented Dickey-Fuller Test

data:  diff(fill_data$move)
Dickey-Fuller = -7.7499, Lag order = 14, p-value = 0.01
alternative hypothesis: stationary
```

p値が大変小さくなりました。
つまり、差分系列は「単位根を持つ」という仮説が棄却され、定常である可能性が少しだけ見えてきました。




#周期は・・・見つかりませんでした

気持ちよく周期が見つけられれば、予測モデルを作成する際にありがたいのです。
周期を間違えるとめちゃくちゃな結果になりますので・・・

#次回は予測モデルの作成について