Tidyverts による時系列予測モデリング
Tags: azure, microsoft, forecasting, machinelearning
Published: 2020-04-21 11:12:06

R を使って、基本的な時系列予測モデルの構築と予測値算出を行っていきます。本投稿では、R の時系列モデリングのライブラリ `tidyverts` を中心に利用します。

### tidyverts とは？
[tidyverts](https://tidyverts.org/) は、大きく下記の3つのライブラリから構成されます。<br/>
[tsibble](https://tsibble.tidyverts.org/) : 時系列データに特化たデータ型を提供

<img align="left" src="https://tidyverts.org/images/tsibble.png">
[fable](https://fable.tidyverts.org/) : 時系列予測のモデルを提供

<img align="left" src="https://tidyverts.org/images/fable.png">
[feasts](https://feasts.tidyverts.org/) : 時系列データの統計処理や特徴量抽出の機能を提供

<img align="left" src="https://tidyverts.org/images/feasts.png">


### Microsoft Forecasting Best Practice

Microsoft が Forecasting の Best Practice 集を公開しています。
{% github microsoft/forecasting %}

このリポジトリでも Tidyverts が利用されていますので、是非ご確認ください。

需要予測などのテーマでは、例えば "店舗" ✖️ "ブランド" ごとにモデルを作る必要があったりします。そのために並列処理の実装が欠かせません。本投稿では言及しませんが、このリポジトリでは、Python では `Ray` 、R では `{parallel}` を利用して並列処理を実装しています。



### 利用するデータ
Prophet のチュートリアルで使用されている Peyton Manning の Wikipedia の閲覧者数の時系列データを利用します。下記 URL から CSV データをダウンロードしておきます。

[example_wp_log_peyton_manning.csv](https://github.com/facebook/prophet/blob/master/examples/example_wp_log_peyton_manning.csv)


### 分析環境

R Studio や Visual Studio Code で構いません。私は Azure Machine Learning で提供している PaaS の R Studio を利用しています。

[Azure Machine Learning](https://azure.microsoft.com/ja-jp/services/machine-learning/)

[R Support in Azure Machine Learning | AI Show | Channel9](https://channel9.msdn.com/Shows/AI-Show/R-in-Azure-Machine-Learning)
{% youtube 9_alql02cT0 %}




 
### ライブラリのロード

利用するライブラリをロードします。

```{r}
library(fable)
library(fable.prophet)
library(tsibble)
library(tsibbledata)
library(feasts)
library(lubridate)
library(dplyr)
library(ggplot2)
```


### データ準備と探索

CSV ファイルのパスを指定して、データをロードします。

```{r}
df <- read.csv("example_wp_log_peyton_manning.csv", colClasses=c("Date","numeric")) 
```

df に格納されているデータを `tsibble`形式に変換し可視化します。"%>%" は R 言語で利用されるパイプ処理を意味しています。autoplot を適用すると可視化まで簡単にできます。

```{r}
df %>% as_tsibble %>% autoplot
```
![](https://user-images.githubusercontent.com/3815677/79829639-9096c600-83de-11ea-8d09-342ab0764ff0.png)
```{r}
df %>% as_tsibble %>% fill_gaps
```
[`fill_gaps`](https://tsibble.tidyverts.org/reference/fill_gaps.html) を利用して、データに欠損している日時項目を補完します。また `tidyr` の [`fill`](https://tidyr.tidyverse.org/reference/fill.html) で欠損値を補完します。１つ先の時系列データ値を用いて補完する場合は、`.direction  = down` と指定します。

```{r}
df_filled <- df %>% as_tsibble %>% fill_gaps %>% tidyr::fill(y, .direction = "down")
```
次に `feasts` ライブラリの [`STL`](https://feasts.tidyverts.org/reference/STL.html) を用いて、時系列データを "季節" "トレンド" "不規則成分" に分解します。

```{r}
dcmp <- df_filled %>%
    model(STL(y ~ season(window = Inf ))) %>% components
dcmp
```
可視化して特徴を確認します。Weekly の季節性トレンドは細かすぎて黒塗りに見えちゃっているようです。
```{r}
dcmp %>% autoplot
```
![](https://user-images.githubusercontent.com/3815677/79829643-955b7a00-83de-11ea-9ea3-c6d2aa2a66b7.png)

### 時系列モデルの構築と予測値算出

まずは非常にシンプルなモデルを試します。データには季節性などが確認されるものの、これらのモデルは考慮してくれないことが分かります。

- MEAN : 平均値
- NAIVE(=RW): ランダムウォークモデル


ここでは、ドリフトを考慮するモデルも追加します。このデータは後半にトレンドが減少している傾向にあるので、このモデルも減少方向に予測していることが確認できます。
```{r}
df_filled %>% 
    model(
        mean=MEAN(y),
        naive=NAIVE(y),
        drift=RW(y ~ drift()),  # ドリフトを考慮
        .safely=FALSE) %>% 
    forecast(h="2 years") %>% 
    autoplot(filter(as_tsibble(df_filled),year(ds) > 2007), level = NULL)
```
![](https://user-images.githubusercontent.com/3815677/79829661-9ee4e200-83de-11ea-93c7-63a2c8063f33.png)

次に季節性を考慮します。

- SNAIVE : 季節性を考慮するランダムウォークモデル

1年周期の季節性があるデータなので、過去 1 年前のデータを予測とするモデルを採用します。

```{r}
df_filled %>% 
    model(
        mean=MEAN(y),
        naive=NAIVE(y),
        drift=RW(y ~ drift()),  # ドリフトを考慮
        snaive = SNAIVE(y ~ lag("year")), # 季節性を考慮
        .safely=FALSE) %>% 
    forecast(h="2 years") %>% 
    autoplot(filter(as_tsibble(df_filled),year(ds) > 2007), level = NULL)
```
![](https://user-images.githubusercontent.com/3815677/79829667-a1dfd280-83de-11ea-8776-df55021c95a4.png)

次に指数平滑化を使用します。

- ETS : 指数平滑化

1年周期を適用してみます。
```{r}
df_filled %>% 
    model(
        est=ETS(y ~ season("A", period= "1 year")),
        .safely=FALSE) %>% 
    forecast(h="1 year") %>% 
    autoplot(filter(as_tsibble(df_filled),year(ds) > 2015), level = NULL)
```
すると、エラーが出てきました。

> Error: Seasonal periods (`period`) of length greather than 24 are not supported by ETS.

最大 24 サイクルしか考慮してくれないようです。このデータは Daily のデータなので、年周期の季節性の考慮は難しそうです。

週周期の季節性があると想定して、モデル構築と予測値算出を行います。

```{r}
df_filled %>% 
    model(
        est=ETS(y ~ season("A", period= "7 days")),
        .safely=FALSE) %>% 
    forecast(h="30 days") %>% 
    autoplot(filter(as_tsibble(df_filled),year(ds) > 2014), level = NULL)
```

![](https://user-images.githubusercontent.com/3815677/79829671-a5735980-83de-11ea-8d9a-a97b7291f9ad.png)

次に ARIMA モデルを利用します。ETS ではうまく動きませんでしたが、ARIMA では、`period` や `D=1` を設定することで年周期の季節性も考慮してくれます。

- (S)ARIMA : (季節)自己回帰和分移動平均

```{r}
df_filled %>% 
    model(
        arima=ARIMA(y ~ pdq(1,1,3) + PDQ(0,1,0,period="365 days")),
        .safely=FALSE) %>%
        forecast(h="2 years") %>% 
    autoplot(filter(as_tsibble(df_filled),year(ds) > 2007), level = NULL)
```
![](https://user-images.githubusercontent.com/3815677/79851597-81c10b00-8400-11ea-8d36-05696f099989.png)

次に、Facebook の Prophet を利用します。`fable.prophet` というライブラリが公開されており、いままで使ってきた `fable` でも Prophet のモデルを組み込めるようになっています。詳細は [fable.prophet](https://github.com/mitchelloharawild/fable.prophet) からご確認ください。

```{r}
df_filled %>% 
    model(
        prophet(y),
        .safely=FALSE) %>% 
    forecast(h="2 years") %>% 
    autoplot(filter(as_tsibble(df_filled),year(ds) > 2007), level = NULL)
```
![](https://user-images.githubusercontent.com/3815677/79836756-8af3ad00-83eb-11ea-8de9-8536f9bd7d41.png)


以上です。