Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
52
Help us understand the problem. What is going on with this article?
@mshinoda88

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

More than 1 year has passed since last update.

はじめに

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

参考資料

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
(出典: 未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge)

ARモデル(自己回帰モデル, autoregressive model)

ARモデルは図のように時間の変化に対し規則的に値が変化する、最も単純な時系列モデル。
ar.jpg
直前のp個の値と相関のあるモデルをAR(p)と表現します。
ARモデルはある時点のデータがそれ以前のデータで回帰的に推定できるモデルです。

MAモデル(移動平均モデル、Moving Average model)

MAモデルは図のように時間の変化に対し不規則に値が変化するけど、
ある区間での変動が一定であるようにモデルを考えます。
ma.jpg
直前のq個の値の誤差の影響を受けるモデルをMA(q)と表現します。
MAモデルは過去の誤差に影響されるモデルです。

ARMAモデル(自己回帰移動平均モデル)

ARモデルとMAモデルの2つ組み合わせて定式化したモデル。
ARモデルに関連する次数をpとMAモデルに関連する次数をqとすると、ARMA(p,q)と表すことができます。

定常性

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

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

SARIMAモデル(季節自己回帰和分移動平均モデル)

sarima.jpg
(出典: 未来を予測するビッグデータの解析手法と「SARIMAモデル」 - DeepAge)

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

データを準備する


# 基本のライブラリを読み込む
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()

スクリーンショット 2017-08-27 23.37.35.png

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

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


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

Dickey-Fuller test

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


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

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)

エラーが出てしまいました。
標準偏差は大きく変動していないけど、移動平均は緩やかに下落している模様。

対数をとって傾向を確認

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


ts_log = np.log(ts)
plt.plot(ts_log, label='rating(log)')
plt.title('rating time series data')
plt.legend(loc='best')

ダウンロード (2).png

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


# 傾向(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

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

SARIMAモデルの最適なパラメータを探す

SARIMA モデルを以下のライブラリでオブジェクト化してみるが、
http://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html


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は年間を四半期に分けているので分割数4と指定しています。


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モデルを作成する。


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.

# 残差のチェック
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

最後に予測

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


# 予測
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

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

今回の解析notebook

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

52
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
mshinoda88
公共系インフラエンジニア、フロントエンジニアを経て、医療系クラウドサービスを展開。医療系システムコンサルタント等を経てAIベンチャーにてデータサイエンティスト。大手から中小まで数多くのシステム開発プロジェクトで開発統括、プロマネを経験。 基盤設計、統計学、機械学習、深層学習、組織論、リーダーシップ論、心理学、事業戦略論等をテーマに、社内外で講師の経験多数。

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
52
Help us understand the problem. What is going on with this article?