ソシャゲ(スクフェス)のボーダーを統計的に予測する方法-非定常時系列データ分析-

******************************
追記
2014/06/20
・少しだけ加筆
2014/06/19
・Rコマンドの部分,一部文字化けしていたので修正
・少しだけ加筆
******************************

 

という訳で,久しぶりに小ネタ.

金融データを用いても良いんだけれど,
というか,その方がネタとしては王道だと思うんだけど,
あえてソシャゲ等の緩いものを本気で考える方が面白いなと思ったのと,

・ イベント期間が決まっている
・ データを入手しやすい
・ 時系列データの傾向が明らか(基本的に右肩上がりのトレンドを持つ)
・ 周期性を予測しやすい(1日の中でアクティブ率が高い時間,低い時間等)

等々,時系列データ解析において最も大事な「仮定(仮説)」を立てやすいので採用.
もう一つは,ファイナンスデータについては色々と持論を持っていまして,
主に経済物理学的観点から,この手の解析は本質的には無意味だと考えているので.
(最も金融データの様な非定常時系列データについても同じ方法論で考える事はできるけど,
仮にそれで有意味な結果が得られるのであれば,僕は今頃大金持ちです\(^o^)/)

●最初に注意事項
・ あくまでネタ記事なので,「非定常時系列データ分析」を真剣に学びたい,やりたい方は参考文献を参照して下さい
(無責任ではありますが,割と適当に書いていく予定なので一々突っ込まれても困る..と先に言い訳)
・ 最小努力の法則がモットー
(省ける処理は省く;一々検定しない,仮定できそうなら勝手に仮定していく)
(トレンドの有無を調べるために重回帰分析を行って回帰係数の有意性をみる,
データが単位根過程であるかどうかADF検定をする等々はせずに,
グラフみてトレンドありそうならトレンドがあると考えるし,
周期があったら周期仮定するし,差分とって弱定常性仮定できそうならしちゃう
(単位根過程とみなしちゃう)
最も今回の様な1次時系列であれば実際的に大きな問題は無いと思う)
・ 処理は全てRで行っていますがコマンドは載せたり載せなかったり適当で
(詳細・・というかちゃんとしたやり方は参考文献を参照して頂ければ)

解析対象のデータは以下から拝借した.
また,何回か欠損値があったので(2014/6/9 16:35とか2014/6/13 00:35とか),
データの欠損部分は適当に補間した.

ラブライブ! スクールアイドルフェスティバルイベントポイント掲載所二代目

次のイベントが始まったら消えちゃうので,
以下は今回使用したデータのスクリーンショット(表は省略).

20130607_累積値

20140607_時速

さて,「予測」について話をする訳なので,イベント開始からイベント終了までの
全データを眺めていても,予測も糞も無いじゃんって事で仕方無いんだけど,
予測をする際に最も大事なのは,データの傾向,性質を掴む事.
なので,とりあえずはせっかくグラフがある事なので,グラフを眺めてみる.
(この知見は勿論今後の解析に役立つし,おおよそこの手のデータでは同じ事が云える)
時系列データの解析において,最も大事な事は,
そのデータのボックスプロットやトレンドグラフをみる事.
人に説明したり,纏める際には数値として見せられる事が望ましい(説得力が増す)ので,
統計学的検定を要するけども,実務的には,検定よりもデータを
視覚的に捉える方が遥かに大事だと思う.

逆に,検定で満足してしまうと,その方が過誤に陥るリスクが高いし,
十分に理解して検定するので無ければ,寧ろしない方がマシだとさえ思う.
αエラーやβエラーの意味p値の意味が分からずにやる位ならしない方がマシ)

対して,人間の脳みそというのは偉大なもので,
グラフをみるだけで,数理化が困難な事でも分かってしまうからね.
特にファイナンスデータを数理的に理解しようとすればする程,
人の直観(行動ファイナンスで云う所のシステム1)がいかに偉大かって感じる.
勿論,それは誤りに繋がる事も往々にしてあるけど..

それでも,オリジナルデータをまずは大切にして,それを眺める事がいかに大切かと.

話が逸れたけども,実際に累積値のグラフをみると,
明らかに右肩上がりのトレンドがある事が分かる.
(実際にトレンドの有無を判断するには重回帰分析を行って回帰係数の有意性をみるが省略)
つまり,この時系列データは非定常データであり,このままだと解析に向かない.
非定常データとは定常性を仮定できないデータを指し,
定常性とは平均,自己共分散が時間に依存しないという性質の事.
簡単に言えば,もし非定常であれば,時系列データは一度しか観測できないにも関わらず,
各統計量が時間に依存するという事になれば,期待値を求めるのは困難だし,
分散や自己共分散を求めるのは不可能になってしまう.
(何かしらの強い仮定(時間構造をモデル化する等)を置く必要が出てくる)
また,トレンドを持つデータ=平均回帰性が無い訳だから,中心極限定理が成立しない.

次にイベント時速グラフをみてみる.
このデータは1時間おきのデータになっているので,
時速グラフというのは,要は差分データのグラフを表している.
これをみると,例えば赤線をみてみると,おおよそ(みた感じ)平均が200で,
そこから上下変動している様にみえる(イベント開始直後と終了間際は一旦おいておく).

つまり,時間に依存せず平均と共分散を算出できそうにみえる為,
この時系列データは単位根過程であり,差分を取れば定常過程とみなせそうだと考えられる.
(ここで本来ならADF検定を行って単位根過程がどうか調べるが省略)

また,イベント初日とイベント最終日では例外的に大きな変動がある事,
差分の動意について,およそ1日の周期を持っていそう,という事が分かる.

簡単に,累積値グラフと時速グラフの自己相関を求めてみると,

acf_x1

acf_diff_x1

累積値はその名の通り値が累積していく為,ラグを取っても以前の値が含まれるので自己相関構造を持つ.
その為,自己相関を取るとこの様なグラフになる.イコールトレンドのあるグラフとも云える.
そして,時速グラフの自己相関グラフをみると,ラグ0とラグ23近辺で自己相関が高い事が分かる.
この2つのグラフからも,差分を取ると良さ気だという事,24時間の周期性が推測される.

一応,Rでは簡単にトレンドと周期変動を分離する事ができるので,
周期性については後でモデル化する時に大事なので,調べておく.

> xts <- ts(as.numeric(x), start=c(05, 16), frequency=24)
> xts.stl <- stl(xts, s.window="periodic")
> plot(xts.stl)

stl_normal

data:オリジナルデータ(=seasonal+trend+remainder)
seasonal:周期変動成分
trend:トレンド成分
remainder:残差

(自分用memo:
Rには似た様な関数で”decompose()”というものが用意されている.
“decompose()”は移動平均による古典的な手法,
“stl()”はloess を用いて時系列データを観測値=トレンド+周期変動+残差に分解する.
一般的には”stl()”がリーズナブルチョイスか.

この結果からも,明確にupward biasなトレンドがあり,
また,明確に周期変動が存在し,その周期は1日周期である事が分かる.
これは,1日を通してアクティブ率が高い時間帯を考えると,
1日の周期を持ちそうだというのは納得いく話で,確からしい結果だろう.
(単純に計算した為,seasonalが負の値を取っているが,頭の中ではイメージとして
トレンドから-400して全体を+400してやれば,整合性は取れるだろう)

ちなみに,周期変動でピークを取っている時間帯は,およそ00時頃と13時頃の様だ.
日中については残差が大きい為,この結果を以って一概には云えないが(というか書いていて気付いたので後述),
少なくとも,仕事が終わってから,帰宅してからという方が集中する時間帯が21~25時頃で,
その中でもピークタイムが00時頃なんだろう.
(どの時間に接続し難いか,を考えると,これについては割と妥当な結果じゃないかと思う)

13時頃のピークについては,最初は適当に
「JKが休み時間にやりまくってるのかなー」なんて思ったけど,
よくよく考えると,イベント終わりが15:00という事で,
正午~15:00の間に物凄い伸びがあるので,それに引っ張られている可能性が高い.
(ここだけみて「コイツ急にJKとか何言ってるんだ」と思われるの嫌だから補足しておくと,
スクフェスしているラブライバーはJKしかいないってネタがあってだな・・・)

そんな訳で,ちょっと最終日を除いたデータで再度検討してみると,

> x2 <- read.csv("data03.csv", header=F)
> x2 <- x2^T
> x2ts <- ts(as.numeric(x2), start=c(05, 16), frequency=24)
> x2ts.stl <- stl(x2ts, s.window="periodic")
> plot(x2ts.stl)

stl_Trim15

綺麗に日中のお山は消えましたとさ.

Seasonalの値を細かくみると,21:00~03:00に掛けて変動が大きく,
特に23:00~01:00(00:00頃にピーク)に動意が大きい事が分かる.
(意外に仕事帰りのJKが大勢いるという事が推察される)
(また,seasonalが最も小さい値を取るのは15:00頃だった)

さて,これだけでも色々とあーだこーだと議論ができるけど,
今やりたいのはボーダー予測.手元にある時系列データから未来の予測をしたい訳だ.
このデータに関して云えば,最終日を除けば,
強いトレンドと綺麗な周期変動がみられる上,
残差はほぼ200以下という事で,予測は容易い訳だけども,
そもそも知りたいのは最終日,最後の値(ボーダー)であって,
最終日を除けば簡単にモデル化できる,なんて何の意味もない.

それに,必ずしも毎回この様に強いトレンドが得られるとは限らないので,
もう少し一般化して,値を予測する事を考えていきたい.

という訳で,ずらずらと書いてきたので,改めてここまでで分かっている事をまとめると,

・ 強いupward biasなトレンドを持つ→非定常時系列データ
・ このデータは単位根過程に従う様にみえる→差分を取れば定常データとして扱える
・ 綺麗な周期性を持っており,その周期はおよそ1日→周期変動成分を考慮

以上から,ARIMA(自己回帰和分移動平均)モデルを用いて,値を予測しようと思う.
(ARIMA(p, d, q)過程はd階差分系列が定常かつ反転可能なARMA(p, q)過程に従う事と同義)

RでARIMAモデルをモデリングするのは非常に簡単で,
例えば12日までのデータでモデリングをして,
その後イベント終了までの予測値を出すと以下の様になる.
(イベント終了日は15日なので3日前のデータで予測してみるという事)

> library(forecast)
> model3_1 <- auto.arima(x3ts, trace=T, stepwise=F, approximation=F)
> summary(model3_1)
Series: x
ARIMA(1,1,1)(1,0,1)[24] with drift
 
Call: auto.arima(x = x3ts, stepwise = F, trace = T, approximation = F)
 
Coefficients:
ar1     ma1    sar1     sma1     drift
0.6136  0.4070  0.8094  -0.5532  230.0311
s.e.  0.1512  0.1712  0.1141   0.1829   36.6458
 
sigma^2 estimated as 5826:  log likelihood = -1004.65
AIC = 2021.29   AICc = 2021.8   BIC = 2040.25
 
In-sample error measures:
ME       RMSE        MAE        MPE       MAPE       MASE
-6.3629445 76.1089030 40.1591227       -Inf        Inf  0.1821988
 

> plot(forecast(model3_1, level=c(95,99.99),h=72), bty="n", pos=5.625)

ARIMA111_101withdrift

“stepwise=F, approximation=F”にしているけども,
計算を省略した部分で最適なモデルを見逃していないか,
近似的な計算をした部分で怪しい事が起こっていないか,
なんて言う事を一々気にするのも馬鹿馬鹿しいので,
計算コストが許す限りはどちらもFalseにしておくと良い.
Traceの結果は余りにも長いので略.

ちなみに,“stepwise=T, approximation=T”とすると以下の様になる.

> model3_2 <- auto.arima(x3ts, trace=T, stepwise=T, approximation=T)
 
ARIMA(2,1,2)(1,0,1)[24] with drift         : 1726.723
ARIMA(0,1,0) with drift         : 2126.107
ARIMA(1,1,0)(1,0,0)[24] with drift         : 1834.124
ARIMA(0,1,1)(0,0,1)[24] with drift         : 2080.345
ARIMA(2,1,2)(0,0,1)[24] with drift         : 1809.101
ARIMA(2,1,2)(2,0,1)[24] with drift         : 1742.432
ARIMA(2,1,2)(1,0,0)[24] with drift         : 1773.149
ARIMA(2,1,2)(1,0,2)[24] with drift         : 1720.314
ARIMA(1,1,2)(1,0,2)[24] with drift         : 1723.317
ARIMA(3,1,2)(1,0,2)[24] with drift         : 1e+20
ARIMA(2,1,1)(1,0,2)[24] with drift         : 1718.989
ARIMA(1,1,0)(1,0,2)[24] with drift         : 1808.363
ARIMA(2,1,1)(1,0,2)[24]                    : 1716.277
ARIMA(2,1,1)(0,0,2)[24]                    : 1828.327
ARIMA(2,1,1)(2,0,2)[24]                    : 1733.281
ARIMA(2,1,1)(1,0,1)[24]                    : 1715.680
ARIMA(2,1,1)                    : 1869.044
ARIMA(1,1,1)(1,0,1)[24]                    : 1737.108
ARIMA(3,1,1)(1,0,1)[24]                    : 1717.793
ARIMA(2,1,0)(1,0,1)[24]                    : 1719.809
ARIMA(2,1,2)(1,0,1)[24]                    : 1717.641
ARIMA(1,1,0)(1,0,1)[24]                    : 1818.289
ARIMA(3,1,2)(1,0,1)[24]                    : 1718.652
ARIMA(2,1,1)(1,0,1)[24] with drift         : 1721.560
ARIMA(2,1,1)(0,0,1)[24]                    : 1834.745
ARIMA(2,1,1)(2,0,1)[24]                    : 1737.790
ARIMA(2,1,1)(1,0,0)[24]                    : 1791.882
 
Best model: ARIMA(2,1,1)(1,0,1)[24]
 
> summary(model3_2)
Series: x3ts
ARIMA(2,1,1)(1,0,1)[24]
 
Call: auto.arima(x = x3ts, stepwise = T, trace = T, approximation = T)
 
Coefficients:
ar1     ar2     ma1    sar1     sma1
0.4715  0.4484  0.6810  0.8408  -0.5772
s.e.  0.4379  0.4198  0.3451  0.0949   0.1631
 
sigma^2 estimated as 5996:  log likelihood = -1008.73
AIC = 2029.46   AICc = 2029.96   BIC = 2048.41
 
In-sample error measures:
ME       RMSE        MAE        MPE       MAPE       MASE
-1.2071298 77.2125951 39.2779923 -0.1439885  0.7569165  0.1782012

ARIMA211_101

数値だけみても中々判断し難いけど,実際にプロットしてみるとどちらのモデルが妥当かよく分かる.

さて,改めてARIMA(1,1,1)(1,0,1)[24] with driftで予測した値と実際の値を比べてみると,

ARIMA111_101withdrift_and_real

黒実線が実測値で,青実線が予測値,オレンジ色のエリアは95%信頼区間で,
黄色の部分は99.99%信頼区間を表している.
(本当はシックスシグマを出したいけど,この関数はCIは4桁までしか入力できない..)

もう少し分かり易くする為,予測値と実測値だけをピックアップし,
その誤差を棒グラフで表したグラフを以下に載せる.

ARIMA111_101withdrift_and_real_and_error

以上から,最終日である15日までは非常に正確に予測できている事が分かる.
逆に言えば,最終日は完全にそれまでとは異なった挙動を示していて,うまく予測ができていない.
ただし,このケースであれば,99.99%信頼区間(黄色エリア)内には収まっている事から,
最終日にオーバーパフォームする事は,ほぼ確実に分かっている訳なので,
この場合であれば,99.99%信頼区間上限を安全圏とみなす事ができるし,
あるいはより安全を取り,シックスシグマやテンシグマを目安とする事が考えられる.

さて,この程度の予測であれば,以上の様に比較的簡単に,単純に求める事ができる.
ただ,ここから更に精度を上げようとすると,途端に難しい話になっていく.
結局,「最終日の誤差」というのは,金融データの予測の難しさ,と同じ様な議論をする事になる.

要は,ファイナンシャルクライシスやバブルの様な動きがまさにこれで,
この手の時系列データというのは往々にして正規分布には従わず,
パレート分布に従う事になる.
その為,正規分布に従うデータであれば95%信頼区間,
ましてや99.99%信頼区間を超える様な値というのは非常に稀にしか起こり得ないのが,
いとも簡単に起こり得てしまう,という事になる.

その上で,このデータをみていくと,最終日の動きはまさにバブルのそれなので,
カオス性をThresholdにおいて,超えた時に別モデルで議論するのが最も妥当なモデルを組めそう.
ただ,計算コストを考えると,余りその方向で真剣に考えようと思わない.
次に考えられる方法としては,一旦誤差が収束し発散する事,
分散が大きく変動する事に着目し,この様な分散不均一性を持つ時系列データに対しては,
GARCHモデルを適用するという方法が考えられる.
一般的なモデルを構築するという意味では,これが最も王道だろうか.
或いは,カーネル多変量的に考えるならば,リーズナブルチョイスはSVRだろう.
SVRは最も完成された手法の一だと思う.本当に.

まあ,しかし,このデータに関して云えば,そこまでする必要もないか.
もう少しデータを集めて,最終日だけベイズ推定的に補正項を導入すれば
十分に実用的な結果が得られる様に思う.

大事なのは安全圏を,そしてできる限りギリギリのラインを知る事.
でも,最悪ボーダー割ってしまう位なら,石を割った方がマシだろう.
と考えると,ARIMAモデル+95%CI,99.99%CI,テンシグマと,
それまでの誤差をみながらベイズ推定的に補正した値とを比較しながら,
判断するのが一番良いんじゃないだろうか.

という訳で,最後に簡単に纏めると,ARIMAモデルというのは
この手のデータに対してリーズナブルチョイスでありながら,
割とうまい具合に予測する事ができるという事.
ただ,そこから精度を上げていこうと思うと,なかなか難しいって事.
色々な条件の中で,様々な手法を試すという方法もあるけど,
現実的には,ベイズ推定的に補正を掛けるのがお手頃だと思う.
(データの特性や時系列的要因が明らかな場合はね)

スクフェスのボーダー予測位なら,2,3回イベントこなせばかなり良い線いくと思うな.
スコアマッチとマラソンでは挙動が異なるのでそれぞれ2,3周分は欲しい所.
(例えばマラソンの場合,面白い位綺麗に2回変曲点があったり)


例えば,今回の分しかデータは持っていないので以下のグラフはあくまでサンプルだけど,
同イベント3周分のデータをグラフ化すると,以下の様になるだろう.

3回分のデータ_サンプル

このままだとイベント毎にトレンドが失われ具合が悪い.
勿論このままでも解析できない事はないけど,ちょっとの工夫で何とかなりそうだ.
つまり,同イベントについては一意なトレンドがあると
仮定してしまえば良い訳で,すると以下の様なグラフで表せる.

3回分のデータ_サンプル_累積

最終日だけ大きな動意があり,それ以外では綺麗な右肩上がりのトレンドを持つ様にみえる.
(まあ・・今の場合は1つのデータ3つ並べているだけなので当たり前なのですが・・・)
これをトレンドと周期変動に分解してみると,以下の様な感じになる.

> n3data_sample <- read.csv("n3data_sample.csv", header=F)
> n3data_sample <- n3data_sample^T
> n3data_sample_ts <- ts(as.numeric(n3data_sample), frequency=24)
> n3data_sample_ts_stl <- stl(n3data_sample_ts, s.window="periodic")
> plot(n3data_sample_ts_stl)

n3data_sample_24stl

どうしても,イベント最終日の動意が別物になっちゃっていて,
予測される周期変動がそれに引きずられて歪んだり残差が増えたりしているが,
上述した1周分のデータのグラフと比較すると,より綺麗に分離できている事が分かる.

ちな,簡単にARIMAモデル求めて,予測値と実測値のグラフ出してみると以下の通り.

n5data_preandact02

“Unable to check for unit roots”ってワーニングが出てしまった.
単位根検定の結果では,このデータは単位根過程に従うと判断して良さそうだけど,
残念ながら,ちょっとちゃんと考えないとうまくは行かない様だ.

まーでも,ある程度N数を増やせば,最終日の別挙動な大きな動意も含めて,
一意なトレンドと周期変動として解釈できるのではないかなー.
実際のデータについても,こんな感じで進めていけば
うまく説明できるモデルが立てられると思う.

まあ,暫くデータ集めて遊んでみるつもり.


どうせなので,次のイベントデータでは次の5つのモデルを比較してみようと思う.
(Rは処理が遅いので途中で挫折する可能性はあるけど・・・)

・ARIMAモデル(上記のお話)
・GARCHモデル(分散不均一性を考慮したモデル)
・SETARモデル(状態変化を前提とした閾値モデル)
・カルマンフィルタ(今回の様なデータでは優位性はなくARIMAと大きな差は出ない気はする)
(計算コストを考えなければパーティクルフィルタ最強なんだけど)
・SVR(計算コストを考えなければ最も優れた手法の一;僕がSVM,SVR大好きっ子というのはある)

参考文献
[1] 沖本竜義:経済・ファイナンスデータの計量時系列分析.朝倉書店,2010/02.
[2] Rで計量時系列分析:AR, MA, ARMA, ARIMAモデル, 予測|銀座で働くデータサイエンティストのブログ
[3] Rで計量時系列分析:単位根過程、見せかけの回帰、共和分、ベクトル誤差修正モデル|銀座で働くデータサイエンティストのブログ
[4] 時系列解析 実践編|海と魚と統計解析
[5] 時系列分析_実践編|Logics of Blue
[6] R の基本パッケージ base, stats 中の時系列オブジェクトの簡易解説|Rの基本パッケージ中の時系列オブジェクト一覧
[7] Rと時系列(2)|[連載]フリーソフトによるデータ解析・マイニング第34回

カテゴリー: 統計学, 日記 パーマリンク