TV朝日の視聴率推移をSARIMAモデルで予測してみる

## はじめに
やりたいことは python3 で時系列分析を行い、予測までしてみること。
特に季節変動を加味したSARIMAモデルを扱ってみたい。
利用するデータとしては、以下の参考資料で、TV朝日の視聴率の時系列分析をしてみる。
きっかけは２０１６年度４Qで大きく利益率を落としているが、視聴率変動がどう影響しているのか、今後の変動はどう予測できるのかを試行した。

### 参考資料
http://www.tv-asahihd.co.jp/contents/ir_setex/
未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge
https://deepage.net/bigdata/2016/10/22/bigdata-analytics.html

## SARIMA モデルの概要
SARIMAモデルは複数の時系列モデルを複合した手法。
以下の図のような関係があります。

![timeline_model.jpg](https://qiita-image-store.s3.amazonaws.com/0/199265/97a37daf-3525-53ef-a114-5d8233e82c18.jpeg)
（出典: 未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge）

### ARモデル（自己回帰モデル, autoregressive model）
ARモデルは図のように時間の変化に対し規則的に値が変化する、最も単純な時系列モデル。
![ar.jpg](https://qiita-image-store.s3.amazonaws.com/0/199265/3a913f26-511b-0c08-2200-34c6189e7205.jpeg)
直前のp個の値と相関のあるモデルをAR(p)と表現します。 
ARモデルはある時点のデータがそれ以前のデータで回帰的に推定できるモデルです。

### MAモデル（移動平均モデル､Moving Average　model）
MAモデルは図のように時間の変化に対し不規則に値が変化するけど、
ある区間での変動が一定であるようにモデルを考えます。
![ma.jpg](https://qiita-image-store.s3.amazonaws.com/0/199265/99a421a3-1e65-bd07-5814-f3679c63e79f.jpeg)
直前のq個の値の誤差の影響を受けるモデルをMA(q)と表現します。
MAモデルは過去の誤差に影響されるモデルです。

### ARMAモデル（自己回帰移動平均モデル）
ARモデルとMAモデルの２つ組み合わせて定式化したモデル。
ARモデルに関連する次数をpとMAモデルに関連する次数をqとすると、ARMA(p,q)と表すことができます。

### 定常性
ARモデル、MAモデル、ARMAモデルはいずれも定常過程を対象とした時系列モデル。
定常過程は、以下の二つを満たすものを言います。
・平均が時間に依存せずに一定である
・異時点間の共分散が時間差のみに依存する
つまり、定常性のある時系列データは「一定の周期で同程度の変動をしている」といえます。

時系列において、平均値が時間的に常に一定であることが保証されているのは稀。
つまり現実に扱う時系列データは非定常過程であることが多い。
実用に際しては、非定常過程の時系列データを対象とした手法を用いる必要があります。
そこで、SARIMAの出番です。

### SARIMAモデル（季節自己回帰和分移動平均モデル）
![sarima.jpg](https://qiita-image-store.s3.amazonaws.com/0/199265/f1782920-786d-8eba-b190-b3beaa671b49.jpeg)
（出典: 未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge）

最後は少しはしょりましたが、本題にGo.

## データを準備する
```python

# 基本のライブラリを読み込む
import numpy as np
import pandas as pd
from scipy import stats
 
# グラフ描画
from matplotlib import pylab as plt
import seaborn as sns
%matplotlib inline
 
 
# グラフを横長にする
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

import matplotlib as mpl
mpl.rcParams['font.family'] = ['serif']

url = "/work/02_timeseries/tvasahi201701.txt"
file2 = "/work/02_timeseries/tvasahi201701_2.txt"
dataNormal = pd.read_csv(url)
dataNormal.head()
```
<img width="824" alt="スクリーンショット 2017-08-27 23.37.35.png" src="https://qiita-image-store.s3.amazonaws.com/0/199265/45910f59-2f17-a9f0-d5f6-732dd49293d5.png">

全日、ゴールデン、プライム、プライム2の時間帯で視聴率データが揃えられる。
今回はプライムタイムのみのデータを利用して分析してみる。なぜこのデータをチョイスしたかというとトレンドがはっきり現れていて、データ分析、予測の学習にはもってこいだから。

まず、このデータをプロットしてみる。

```python

ts = ((dataNormal.T)[2])

dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m-01')
ts = pd.read_csv(file2, index_col="time", date_parser=dateparse, dtype='float')
plt.plot(ts, label='rating')
plt.title('Viewer rating series data')
plt.legend(loc='best')
```

![ダウンロード.png](https://qiita-image-store.s3.amazonaws.com/0/199265/585b6180-743c-eef0-231a-e964fe8c9666.png)

## Dickey-Fuller test
定常性をチェックするための統計的テストの1つです。
ここで、帰無仮説は、ts が非定常であるということ。
テスト結果は、「検定統計量(Test Statistic)」と1,5,10%の信頼水準の「臨界値(Critical Value)」から構成されます。
「検定統計量」が「臨界値」よりも小さい場合は、帰無仮説を棄却して系列が定常状態と判定できます。


```python

from statsmodels.tsa.stattools import adfuller

# Dickey-Fuller test 結果と標準偏差、平均のプロット
def test_stationarity(timeseries, window_size=12):
    # Determing rolling statistics
    rolmean = timeseries.rolling(window=window_size,center=False).mean()
    rolstd = timeseries.rolling(window=window_size,center=False).std()

    # Plot rolling statistics:
    orig = plt.plot(timeseries, color='blue',label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    # Perform Dickey-Fuller test:
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value',
                                             '#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value (%s)'%key] = value
    print(dfoutput)
```

まずは上記関数を呼び出してみる。
```python

test_stationarity(ts, window_size=4)
```
![ダウンロード (1).png](https://qiita-image-store.s3.amazonaws.com/0/199265/e1a2e4b5-7fb7-6149-c6ed-eb4627602cbc.png)

```
Results of Dickey-Fuller Test:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-8-552a0eb0606f> in <module>()
----> 1 test_stationarity(ts, window_size=4)

<ipython-input-7-e5ce27c1c7e1> in test_stationarity(timeseries, window_size)
     17     # Perform Dickey-Fuller test:
     18     print('Results of Dickey-Fuller Test:')
---> 19     dftest = adfuller(timeseries, autolag='AIC')
     20     dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value',
     21                                              '#Lags Used','Number of Observations Used'])

/root/anaconda3/lib/python3.6/site-packages/statsmodels/tsa/stattools.py in adfuller(x, maxlag, regression, autolag, store, regresults)
    218 
    219     xdiff = np.diff(x)
--> 220     xdall = lagmat(xdiff[:, None], maxlag, trim='both', original='in')
    221     nobs = xdall.shape[0]  # pylint: disable=E1103
    222 

/root/anaconda3/lib/python3.6/site-packages/statsmodels/tsa/tsatools.py in lagmat(x, maxlag, trim, original, use_pandas)
    373     if xa.ndim == 1:
    374         xa = xa[:, None]
--> 375     nobs, nvar = xa.shape
    376     if original in ['ex', 'sep']:
    377         dropidx = nvar

ValueError: too many values to unpack (expected 2)
```
エラーが出てしまいました。
標準偏差は大きく変動していないけど、移動平均は緩やかに下落している模様。

## 対数をとって傾向を確認
次に対数をとって傾向を確認して見ましょう。

```python

ts_log = np.log(ts)
plt.plot(ts_log, label='rating(log)')
plt.title('rating time series data')
plt.legend(loc='best')
```
![ダウンロード (2).png](https://qiita-image-store.s3.amazonaws.com/0/199265/e3ae930e-a272-d220-b63a-e6564f7e2b7e.png)

次に対数データについて、傾向、季節性、残差を眺めてみる。

```python

# 傾向(trend)、季節性(seasonal)、残差(residual)に分解してモデル化する。
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log)

trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid


# オリジナルの時系列データプロット
plt.subplot(411)
plt.plot(ts_log, label='Original')
plt.legend(loc='best')
 
# trend のプロット
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
 
# seasonal のプロット
plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
 
# residual のプロット
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()
```
![ダウンロード (3).png](https://qiita-image-store.s3.amazonaws.com/0/199265/004ebcd5-ff9d-0007-a3d6-b6f99c15e55a.png)

秋にかけて視聴率が下がり、冬に視聴率が上がるというのが毎年の季節性のある傾向。
トレンドとしては緩やかに下がってきているのが見て取れます。

## SARIMAモデルの最適なパラメータを探す
SARIMA モデルを以下のライブラリでオブジェクト化してみるが、
http://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html

```python

sarima = sm.tsa.SARIMAX(
    ts, order=(p,d,q), 
    seasonal_order=(sp,sd,sq,4), 
    enforce_stationarity = False, 
    enforce_invertibility = False)
```
p,d,q,sp,sd,sq の組み合わせを探してみる。
(sp,sd,sq,4)の最後の４は年間を四半期に分けているので分割数4と指定しています。

```python

import statsmodels.api as sm

# 自動SARIMA選択
num = 0

for p in range(1, max_p + 1):
    for d in range(0, max_d + 1):
        for q in range(0, max_q + 1):
            for sp in range(0, max_sp + 1):
                for sd in range(0, max_sd + 1):
                    for sq in range(0, max_sq + 1):
                        sarima = sm.tsa.SARIMAX(
                            ts, order=(p,d,q), 
                            seasonal_order=(sp,sd,sq,4), 
                            enforce_stationarity = False, 
                            enforce_invertibility = False
                        ).fit()
                        modelSelection.iloc[num]["model"] = "order=(" + str(p) + ","+ str(d) + ","+ str(q) + "), season=("+ str(sp) + ","+ str(sd) + "," + str(sq) + ")"
                        modelSelection.iloc[num]["aic"] = sarima.aic
                        num = num + 1

# モデルごとの結果確認
print(modelSelection)

# AIC最小モデル
print(modelSelection[modelSelection.aic == min(modelSelection.aic)])
```

以下が実行結果。

```
                           model      aic
0    order=(1,0,0), season=(0,0,0)  70.3877
1    order=(1,0,0), season=(0,0,1)  48.7876
2    order=(1,0,0), season=(0,1,0)  40.1175
3    order=(1,0,0), season=(0,1,1)  27.9964
4    order=(1,0,0), season=(1,0,0)   40.728
5    order=(1,0,0), season=(1,0,1)  37.2894
6    order=(1,0,0), season=(1,1,0)  30.3103
7    order=(1,0,0), season=(1,1,1)   29.522
8    order=(1,0,1), season=(0,0,0)  62.5256
9    order=(1,0,1), season=(0,0,1)  43.3847
10   order=(1,0,1), season=(0,1,0)  37.7301
11   order=(1,0,1), season=(0,1,1)  26.0388
12   order=(1,0,1), season=(1,0,0)  41.8186
13   order=(1,0,1), season=(1,0,1)  32.2768
14   order=(1,0,1), season=(1,1,0)  31.8695
15   order=(1,0,1), season=(1,1,1)  28.0252
16   order=(1,0,2), season=(0,0,0)  57.7646
17   order=(1,0,2), season=(0,0,1)  40.7138
18   order=(1,0,2), season=(0,1,0)  36.0716
19   order=(1,0,2), season=(0,1,1)  25.3384
20   order=(1,0,2), season=(1,0,0)  42.2375
21   order=(1,0,2), season=(1,0,1)  28.2212
22   order=(1,0,2), season=(1,1,0)  33.8669
23   order=(1,0,2), season=(1,1,1)  27.4596
24   order=(1,0,3), season=(0,0,0)    57.36
25   order=(1,0,3), season=(0,0,1)  39.7003
26   order=(1,0,3), season=(0,1,0)  37.0603
27   order=(1,0,3), season=(0,1,1)    24.85
28   order=(1,0,3), season=(1,0,0)  45.7587
29   order=(1,0,3), season=(1,0,1)  30.5274
..                             ...      ...
162  order=(3,1,0), season=(0,1,0)  35.6774
163  order=(3,1,0), season=(0,1,1)  26.6998
164  order=(3,1,0), season=(1,0,0)  36.7483
165  order=(3,1,0), season=(1,0,1)  37.3619
166  order=(3,1,0), season=(1,1,0)  27.1641
167  order=(3,1,0), season=(1,1,1)  23.5206
168  order=(3,1,1), season=(0,0,0)  46.0088
169  order=(3,1,1), season=(0,0,1)  42.0283
170  order=(3,1,1), season=(0,1,0)  33.2783
171  order=(3,1,1), season=(0,1,1)  23.3888
172  order=(3,1,1), season=(1,0,0)  33.2987
173  order=(3,1,1), season=(1,0,1)  34.8924
174  order=(3,1,1), season=(1,1,0)  24.7505
175  order=(3,1,1), season=(1,1,1)  25.6202
176  order=(3,1,2), season=(0,0,0)  40.3796
177  order=(3,1,2), season=(0,0,1)  37.4187
178  order=(3,1,2), season=(0,1,0)  35.7897
179  order=(3,1,2), season=(0,1,1)  23.3374
180  order=(3,1,2), season=(1,0,0)  38.7358
181  order=(3,1,2), season=(1,0,1)  37.6919
182  order=(3,1,2), season=(1,1,0)  29.2141
183  order=(3,1,2), season=(1,1,1)  27.3448
184  order=(3,1,3), season=(0,0,0)  41.1869
185  order=(3,1,3), season=(0,0,1)  34.1721
186  order=(3,1,3), season=(0,1,0)  34.8843
187  order=(3,1,3), season=(0,1,1)  23.1797
188  order=(3,1,3), season=(1,0,0)  40.3741
189  order=(3,1,3), season=(1,0,1)  35.8745
190  order=(3,1,3), season=(1,1,0)  23.1581
191  order=(3,1,3), season=(1,1,1)  29.1668

[192 rows x 2 columns]
                            model      aic
59  order=(1,1,3), season=(0,1,1)  20.1616
```
「order=(1,1,3), season=(0,1,1)」が最適なパラメータらしいので、
このパラメータでSARIMAモデルを作成する。

```python

p=1
d=1
q=3
sp=0
sd=1
sq=1

sarima = sm.tsa.SARIMAX(
    ts, order=(p,d,q), 
    seasonal_order=(sp,sd,sq,4), 
    enforce_stationarity = False, 
    enforce_invertibility = False
).fit()
# 結果確認
print(sarima.summary())
```

```
 Statespace Model Results                                
=========================================================================================
Dep. Variable:                              プライム   No. Observations:                   28
Model:             SARIMAX(1, 1, 3)x(0, 1, 1, 4)   Log Likelihood                  -4.081
Date:                           Sat, 26 Aug 2017   AIC                             20.162
Time:                                   10:16:53   BIC                             28.155
Sample:                               04-01-2010   HQIC                            22.605
                                    - 01-01-2017                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.7597   2.03e+06   3.74e-07      1.000   -3.99e+06    3.99e+06
ma.L1         -1.2405      0.000  -5216.846      0.000      -1.241      -1.240
ma.L2         -0.0121      0.003     -4.468      0.000      -0.017      -0.007
ma.L3          0.2525      0.011     22.090      0.000       0.230       0.275
ma.S.L4       -1.0001      0.013    -77.808      0.000      -1.025      -0.975
sigma2         0.0624      0.015      4.218      0.000       0.033       0.091
===================================================================================
Ljung-Box (Q):                        5.47   Jarque-Bera (JB):                 0.47
Prob(Q):                              0.98   Prob(JB):                         0.79
Heteroskedasticity (H):               0.21   Skew:                            -0.43
Prob(H) (two-sided):                  0.11   Kurtosis:                         3.02
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 2.18e+36. Standard errors may be unstable.
```

```python

# 残差のチェック
residSARIMA = sarima.resid
fig = plt.figure(figsize=(12,8))

# 自己相関
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(residSARIMA, lags=4, ax=ax1)

# 偏自己相関
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(residSARIMA, lags=4, ax=ax2)
```
残差の自己相関については、ほぼ問題なくなったことを確認。

![ダウンロード (4).png](https://qiita-image-store.s3.amazonaws.com/0/199265/c325fe2b-7c61-428a-c0d2-5bf654701812.png)

## 最後に予測
長々と色々見てきましたが、SARIMA モデルでどれだけ予測できるのかを見てみます。
元データは２０１０年1Qから2016年4Qですが、この元データと予測データ(2015年1Q-2017年4Q)までを同時プロットしてみます。

```python

# 予測
ts_pred = sarima.predict('2015-01-01', '2017-12-01')

# 実データと予測結果の図示
plt.plot(ts, label='original')
plt.plot(ts_pred, label='predicted', color='red')
plt.legend(loc='best')

print(ts_pred)
```
```
2015-01-01    11.601939
2015-04-01    10.927707
2015-07-01    10.164144
2015-10-01    11.711078
2016-01-01    11.266165
2016-04-01    10.183324
2016-07-01     9.949642
2016-10-01    11.106567
2017-01-01    10.775023
2017-04-01     9.803198
2017-07-01     9.344881
2017-10-01    10.452862
Freq: QS-OCT, dtype: float64
```
![ダウンロード (5).png](https://qiita-image-store.s3.amazonaws.com/0/199265/e3f6e3bf-4ba4-0897-ffbf-a4042becf0a4.png)


２０１５-２０１６年度のデータはまずまずの精度で予測できていることが分かります。
２０１７年度は秋口にかけて大きく下落し、冬に大きく上昇することが予測できそう。
トレンドを加味して予測するのはさほど難しくないですね。
今回はcost function による厳密なcost分析を行っていませんが、
今後のアンサンブル解析時の課題としたいと思います。

## 今回の解析notebook
今回の詳細データ解析過程については notebook を参照してください。
https://github.com/mshinoda88/python/blob/master/20170822-01_tvasahi.ipynb