ソシャゲ（スクフェス）のボーダーを統計的に予測する方法-非定常時系列データ分析-｜粉末@それは風のように (日記)
の続き．
現在進行形のイベントデータを元に，ボーダー予測をしてみる．
データは前回と同じ様に以下から拝借．
欠損値があれば，補間している．
ラブライブ！ スクールアイドルフェスティバルイベントポイント掲載所二代目
今回は主に，色んな手法で予測をしてみようと思う．
なるべく，結果に対しての考察や解釈は控えるつもり．
というより，いろんな手法で解析する事に重きを置き過ぎて，
あーだこーだ考察したり，理屈を説明する余力は全く無かった……．
マラソンイベでデータ取るのは今回が初めてなので，
どの結果が最も妥当か等の考察も行わない（行えない）．
各手法のグラフをみて，イメージを膨らませられれば良いかなと．
後，ポイントは26日（木）にEX曲追加される訳だけど，
その時に確実にそれまでと増分率が変わる事が予想される．
それによって，いかに初期モデルでの予測結果が実測結果と異なるか，とか，
そういう事（インパクト）が何時あるかとか分かっているのに，
いかに最終値を精度良く予測するのが困難かとか，
ましてファイナンスデータなんて・・・とか，
なんかそういう事に思いを馳せると面白いかも．
詳しい理屈，説明は参考文献を参照．ただし，自分用memoとして，
この記事をみただけで同じ事ができる程度には丁寧に書いていきたい（願望）．
（今まで主にMATLABやC++/JAVAでゴニョゴニョやってきたのでRの備忘録的な）
思いつきで手法をチョイスしているので，何故その手法をとか特に理由は無い．
（本当はその辺りの事を考えるのが大切だし，また最も面白い部分なんだけど）
・解析するデータ
－6500位のスコア
－32,500位のスコア
・手法
－線形回帰モデルによる予測（lm()関数）
－一般化線形モデルによる予測（glm()関数）
－一般化線形混合モデルによる予測（library(nlme)，lme()関数）
－非線形回帰モデルによる予測（nls()関数）
－非線形混合回帰モデルによる予測（library(nlme），nlme()関数)
－一般化加法モデル（GAM）による予測（library(mgcv)，gam()関数）
－指数平滑法（ETS）による予測（library(forecast)，ets()関数）
－ARIMA（自己回帰和分移動平均）モデルによる予測（前回記事[1]の話）
－ARFIMA（自己回帰実数和分移動平均）モデルによる予測（library(fracdiff)，fdGPH()やfracdiff()）
－GARCH（一般化自己回帰条件付分散不均一）モデルによる予測（library(fSeries)，garchfit()関数）
－SETAR（自己励起型閾値）モデルによる予測（library(tsDyn)，setar()関数）
－カルマンフィルタによる予測（予定）
－パーティクルフィルタによる予測（予定）
－SVR（サポートベクター回帰）モデルによる予測（予定）
●2014/6/20(金) 16:00～2014/6/25（水） 18:35のデータから予測
・線形回帰モデルによる予測
－6500位のスコア
Rコマンド
xy <- read.csv("data_6500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
xy_lmmodel <- lm(y~x-1, xy)#単回帰（線形回帰）モデル：y=ax（切片は0）
summary(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel, se.fit=T, interval="confidence")
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_lmpred_se <- 2*(sefit$coefficients[2])
new_x <- seq(0, 240,1)
new_y <- new_x*xy_lmmodel_coefficients#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_lmpred_se
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_lmmodel) Call: lm(formula = y ~ x - 1, data = xy) Residuals: Min 1Q Median 3Q Max -600.16 -80.58 251.66 409.46 842.50 Coefficients: Estimate Std. Error t value Pr(>|t|) x 158.813 1.035 153.5 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 384.1 on 74 degrees of freedom Multiple R-squared: 0.9969, Adjusted R-squared: 0.9968 F-statistic: 2.356e+04 on 1 and 74 DF, p-value: < 2.2e-16
> par(mfrow=c(2,2)) > plot(xy_lmmodel)
決定係数，自由度調整済み決定係数から，データへの当て嵌まりは良いが，
残差にトレンドがみられる為，どんどん誤差が大きくなっていく事が容易に想像できる．
また，残差に周期変動が見られる事から，周期変動を含めてモデリングすると良さそう．
残差の正規Q-Qプロットの結果から，残差は正規分布に従うと考えて良さそう．
[235,] NA 37162.2221 36678.0141 37646.4301
[236,] NA 37321.0350 36834.7577 37807.3123
[237,] NA 37479.8479 36991.5014 37968.1945
[238,] NA 37638.6608 37148.2450 38129.0766
[239,] NA 37797.4738 37304.9887 38289.9588
[240,] NA 37956.2867 37461.7323 38450.8410
[241,] NA 38115.0996 37618.4760 38611.7232
予測値：38115.0996±496.6236
－32,500位のスコア
xy <- read.csv("data_32500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
xy_lmmodel <- lm(y~x-1, xy)#単回帰（線形回帰）モデル：y=ax（切片は0）
summary(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel, se.fit=T, interval="confidence")
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_lmpred_se <- 2*(sefit$coefficients[2])
new_x <- seq(0, 240,1)
new_y <- new_x*xy_lmmodel_coefficients#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_lmpred_se
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_lmmodel) Call: lm(formula = y ~ x - 1, data = xy) Residuals: Min 1Q Median 3Q Max -246.001 5.333 64.000 146.375 280.292 Coefficients: Estimate Std. Error t value Pr(>|t|) x 67.9583 0.3834 177.2 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 142.3 on 74 degrees of freedom Multiple R-squared: 0.9977, Adjusted R-squared: 0.9976 F-statistic: 3.142e+04 on 1 and 74 DF, p-value: < 2.2e-16
> par(mfrow=c(2,2)) > plot(xy_lmmodel)
決定係数，自由度調整済み決定係数から，データへの当て嵌まりは良いが，
残差にトレンドがみられる為，どんどん誤差が大きくなっていく事が容易に想像できる．
また，残差に周期変動が見られる事から，周期変動を含めてモデリングすると良さそう．
残差の正規Q-Qプロットの結果から，残差は正規分布に従うと考えて良さそう．
[235,] NA 15902.25290 15584.65350 16219.85231
[236,] NA 15970.21125 15651.37926 16289.04324
[237,] NA 16038.16959 15718.10502 16358.23417
[238,] NA 16106.12794 15784.83078 16427.42510
[239,] NA 16174.08628 15851.55654 16496.61603
[240,] NA 16242.04463 15918.28230 16565.80696
[241,] NA 16310.00297 15985.00806 16634.99789
予測値：16310.00297±324.99492
－一般化線形モデルによる予測（glm()関数）
線形回帰モデルでは残差が正規分布に従うと仮定している．
一般化線形モデルでは，リンク関数により正規分布を拡張した分布族に対応させ，
非線形現象を線形モデルで扱える様にしようと考える方法．
（観測空間を線形分離可能な特徴空間に変換して解くイメージ；そして始まるカーネルトリック）
前述の通り，このデータでは残差が正規分布に従うと仮定できるため，
ここではfamily=gaussian（glm()のデフォルト値）として，
周期変動を含めたモデリングを行う（やっている事はlm()で重回帰と同じ）．
具体的には，日内変動を考慮して，0~23時というファクターをおいて，
目的変数：スコア
説明変数：日時，経過時間
という風に考えてみる．
– 6,500位のスコア
Rコマンド
xy <- read.csv("data_6500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
atemp <- rep(1:24, 4)
a <- atemp[1:length(xy[,1])]
xy <- cbind(xy, a)
xy_glmmodel <- glm(y~x-1+a, xy, family=gaussian)#lm(y~x-1+a, xy)と同じ
summary(xy_glmmodel)
xy_glmpred <- predict(xy_glmmodel, se.fit=T, interval="confidence")
xy_glmmodel_coefficients <- xy_glmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_glmpred_se <- 2*(sefit$coefficients[2])#2SE（95%信頼区間を求める;自由度は考慮しない）
new_x <- seq(0, 240,1)
new_a <- rep(1:24,11)
new_a <- new_a[1:length(new_x)]
new_y <- new_x*xy_glmmodel_coefficients[1] + new_a*xy_glmmodel_coefficients[2]#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_glmpred_se#それに信頼区間を付ける為標準偏差求める
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_glmmodel) Call: glm(formula = y ~ x - 1 + a, family = gaussian, data = xy) Deviance Residuals: Min 1Q Median 3Q Max -526.5 -132.8 136.5 373.6 706.4 Coefficients: Estimate Std. Error t value Pr(>|t|) x 153.697 1.547 99.335 < 2e-16 *** a 19.672 4.735 4.154 8.76e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 120961.1) Null deviance: 3487076739 on 75 degrees of freedom Residual deviance: 8830164 on 73 degrees of freedom AIC: 1094.6 Number of Fisher Scoring iterations: 2
> par(mfrow=c(2,2)) > plot(xy_glmmodel)
[235,] NA 36338.94624 36021.346833 36656.54564
[236,] NA 36512.31592 36193.483934 36831.14792
[237,] NA 36685.68561 36365.621036 37005.75019
[238,] NA 36859.05530 36537.758137 37180.35246
[239,] NA 37032.42499 36709.895238 37354.95473
[240,] NA 37205.79467 36882.032340 37529.55701
[241,] NA 36907.02737 36582.032449 37232.02229
累積値が減るという事はあり得ないので，モデルとしては宜しくない．
– 32,500位のスコア
xy <- read.csv("data_32500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
atemp <- rep(1:24, 4)
a <- atemp[1:length(xy[,1])]
xy <- cbind(xy, a)
xy_glmmodel <- lm(y~x-1+a, xy)
summary(xy_glmmodel)
xy_glmpred <- predict(xy_glmmodel, se.fit=T, interval="confidence")
xy_glmmodel_coefficients <- xy_glmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_glmpred_se <- 2*(sefit$coefficients[2])#2SE（95%信頼区間を求める;自由度は考慮しない）
new_x <- seq(0, 240,1)
new_a <- rep(1:24,11)
new_a <- new_a[1:length(new_x)]
new_y <- new_x*xy_glmmodel_coefficients[1] + new_a*xy_glmmodel_coefficients[2]#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_glmpred_se#それに信頼区間を付ける為標準偏差求める
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 20000), pos=0)
結果：
> summary(xy_glmmodel) Call: lm(formula = y ~ x - 1 + a, data = xy) Residuals: Min 1Q Median 3Q Max -271.85 -42.31 77.98 148.17 245.91 Coefficients: Estimate Std. Error t value Pr(>|t|) x 66.511 0.601 110.665 < 2e-16 *** a 5.564 1.839 3.025 0.00343 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 135.1 on 73 degrees of freedom Multiple R-squared: 0.9979, Adjusted R-squared: 0.9979 F-statistic: 1.744e+04 on 2 and 73 DF, p-value: < 2.2e-16
> par(mfrow=c(2,2)) > plot(xy_glmmodel)
[235,] NA 15669.40529 15546.036756 15792.77383
[236,] NA 15741.48074 15617.633416 15865.32806
[237,] NA 15813.55618 15689.230076 15937.88229
[238,] NA 15885.63163 15760.826736 16010.43653
[239,] NA 15957.70708 15832.423396 16082.99076
[240,] NA 16029.78252 15904.020056 16155.54499
[241,] NA 15968.32318 15842.081924 16094.56443
予測値：15968.32318±126.24125
－一般化線形混合モデルによる予測（library(nlme)，lme()関数）
この辺り，ちょっと理解が浅いのであれなんだけど……この辺は非常に面白い話．
信頼区間の算出方法が分からなかったので，sd(DATA[,2])/sqrt(length(xy[,1]))より求めている．
－6500位のスコア
xy <- read.csv("data_6500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
atemp <- rep(1:24, 4)
a <- atemp[1:length(xy[,1])]
xy <- cbind(xy, a)
gxy <-groupedData(y ~ x | a, data =xy)
xy_lmmodel <- lme(y~x, gxy)
summary(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel)
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
xy_coeff <- xy_lmmodel_coefficients$fixed #[1]が切片，[2]がxの係数
xy_random <- xy_lmmodel_coefficients$random #$a[a,1]が切片，$a[a,2]が係数
new_x <- seq(0, 240,1)
new_a <- rep(1:24,11)
new_a <- new_a[1:length(new_x)]
new_y <- xy_coeff[1] + new_x*xy_coeff[2] + xy_random$a[new_a,1] + new_x*xy_random$a[new_a, 2]#将来の値を予測
new_CI <- cbind(new_y-2*(sd(DATA[,2])/sqrt(length(xy[,2]))), new_y+2*(sd(DATA[,2])/sqrt(length(xy[,2]))))#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_lmmodel)
Linear mixed-effects model fit by REML
Data: gxy
AIC BIC logLik
941.9017 955.6444 -464.9508
Random effects:
Formula: ~x | a
Structure: General positive-definite
StdDev Corr
(Intercept) 185.414370 (Intr)
x 1.292521 -0.325
Residual 77.896682
Fixed effects: y ~ x
Value Std.Error DF t-value p-value
(Intercept) 625.2525 42.16241 50 14.82962 0
x 146.6253 0.51066 50 287.13171 0
Correlation:
(Intr)
x -0.479
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.94565829 -0.21940089 0.09860771 0.36847152 2.80529837
Number of Observations: 75
Number of Groups: 24
> plot(xy_lmmodel)
23 NA 34926.8760 32568.33750 37285.415
24 NA 35176.9562 32818.41769 37535.495
1 NA 35690.6921 33332.15361 38049.231
2 NA 35350.5221 32991.98357 37709.061
3 NA 35308.8783 32950.33981 37667.417
4 NA 35870.7074 33512.16895 38229.246
予測値：35870.7074±2358.539
– 32,500位のスコア
xy <- read.csv("data_32500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
atemp <- rep(1:24, 4)
a <- atemp[1:length(xy[,1])]
xy <- cbind(xy, a)
gxy <-groupedData(y ~ x | a, data =xy)
xy_lmmodel <- lme(y~x, gxy)
summary(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel)
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
xy_coeff <- xy_lmmodel_coefficients$fixed #[1]が切片，[2]がxの係数
xy_random <- xy_lmmodel_coefficients$random #$a[a,1]が切片，$a[a,2]が係数
new_x <- seq(0, 240,1)
new_a <- rep(1:24,11)
new_a <- new_a[1:length(new_x)]
new_y <- xy_coeff[1] + new_x*xy_coeff[2] + xy_random$a[new_a,1] + new_x*xy_random$a[new_a, 2]#将来の値を予測
new_CI <- cbind(new_y-2*(sd(DATA[,2])/sqrt(length(xy[,2]))), new_y+2*(sd(DATA[,2])/sqrt(length(xy[,2]))))#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_lmmodel)
Linear mixed-effects model fit by REML
Data: gxy
AIC BIC logLik
847.3338 861.0765 -417.6669
Random effects:
Formula: ~x | a
Structure: General positive-definite
StdDev Corr
(Intercept) 42.6747994 (Intr)
x 0.8408221 1
Residual 47.2741676
Fixed effects: y ~ x
Value Std.Error DF t-value p-value
(Intercept) 203.4106 14.020375 50 14.50822 0
x 64.1411 0.311893 50 205.65071 0
Correlation:
(Intr)
x -0.225
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.24126730 -0.36398803 0.03079293 0.54189627 1.97695033
Number of Observations: 75
Number of Groups: 24
>
> plot(xy_lmmodel)
23 NA 15076.3142 12717.77566 17434.853
24 NA 15123.7783 12765.23977 17482.317
1 NA 15120.3367 12761.79817 17478.875
2 NA 15287.8043 12929.26576 17646.343
3 NA 15390.7859 13032.24743 17749.324
4 NA 15593.8971 13235.35857 17952.436
予測値：15593.8971±2358.5389
－非線形回帰モデルによる予測（nls()関数）
Rでは非線形回帰モデルを設計する為の関数として，nls()関数がある．
これは，関数式を自分で指定して与える事ができる．
また，Rに用意された関数式をリンクして非線形回帰モデルを立てる事もできる．
多項式関数poly，ロジスティック関数SSlogis，ワイブル関数SSweibull，勾配関数SSfol等．
ここでは，多項式関数polyを用いて，3次多項式回帰モデルを求める．
・・が，係数の抜き出し方が分からなかったので，多項式求めた後はExcelで図示している．．
－6,500位のスコア
xy <- read.csv("data_6500_2.csv")
xy_nlsmodel <- nls(y~b*x+c*x^2+d*x^3, xy, start=c(b=1, c=1, d=1), trace=T)
summary(xy_nlsmodel)
結果：
> summary(xy_nlsmodel)
Formula: y ~ b * x + c * x^2 + d * x^3
Parameters:
Estimate Std. Error t value Pr(>|t|)
b 209.232852 6.404743 32.668 < 2e-16 ***
c -1.581517 0.266409 -5.936 9.40e-08 ***
d 0.011395 0.002645 4.308 5.12e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 237.7 on 72 degrees of freedom
Number of iterations to convergence: 1
Achieved convergence tolerance: 4.728e-08
>plot(xy_nlsmodel)
黒実線が実測値，赤実線がRによる3次多項式，灰色実線はExcelによる3次多項式．
単に多項式近似を求めるだけなら，Excelの方が楽なんだけど……割と信用できない（後述）．
経験的に考えて，あり得ない値なので，検討するまでもない（経験って大事……）．
－32,500位のスコア
xy <- read.csv("data_32500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
xy_nlsmodel <- nls(y~b*x+c*x^2+d*x^3, xy, start=c(b=1, c=1, d=1), trace=T)#単回帰（線形回帰）モデル：y=ax（切片は0）
summary(xy_nlsmodel)
結果：
> summary(xy_nlsmodel) Formula: y ~ b * x + c * x^2 + d * x^3 Parameters: Estimate Std. Error t value Pr(>|t|) b 79.5368732 2.5899514 30.710 <2e-16 *** c -0.2633746 0.1077307 -2.445 0.0169 * d 0.0009421 0.0010695 0.881 0.3813 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 96.13 on 72 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 2.861e-08
赤色がR，灰色がExcel．Excelの多項式近似は往々にして結果が信用できない．
Excelの統計ツールも怪しい物が多いし（標準のもの；有料パッケージは別），
内部処理がみえないので，検証が困難だし全く信用できない．
Excelでの解析はやらない方が良い（ましてや論文でとか論文そのものが信用できないレベル）．
（Excelで信号処理，統計解析を行っている論文は結果も内容も著者も疑った方が良い）
（ただし有料パッケージは別；割と良いプロダクトはあったりする）
グラフをプロットする位の用途でしか使わない方が良い．
予測値：16830.83
－非線形混合回帰モデルによる予測（library(nlme），nlme()関数)
略
－一般化加法モデル（GAM）による予測（library(mgcv)，gam()関数）
－6,500位のスコア
library(mgcv)
xy <- read.csv("data_6500_2.csv")
gam.model <- gam(y~s(x), xy)
summary(gam.model)
plot(gam.model)
結果：
> summary(gam.model) Family: gaussian Link function: identity Formula: y ~ s(x) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6043.16 9.47 638.2 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(x) 8.838 8.992 12362 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.999 Deviance explained = 99.9% GCV = 7741.1 Scale est. = 6725.6 n = 75
R-sq.(adj) = 0.999 Deviance explained = 99.9%と非常に当て嵌まりが良い．
およそx=36で6,000の上昇と見積もると，最終値は40166.67となる．
（正直，どうしたら良いか分からん・・・）
予測値：40166.67
－32,500位のスコア
xy <- read.csv("data_32500_2.csv")
gam.model <- gam(y~s(x), data=xy)
summary(gam.model)
plot(gam.model)
predict(gam.model, se.fit=T, interval="confidence")
結果：
> summary(gam.model) Family: gaussian Link function: identity Formula: y ~ s(x) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2571.21 3.49 736.8 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(x) 8.88 8.996 17298 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 1 Deviance explained = 100% GCV = 1051.9 Scale est. = 913.34 n = 75
予測値：16736
－指数平滑法（ETS）による予測（library(forecast)，ets()関数）
指数平滑法によるモデリングはARIMAモデルの特殊なケースとして理解できる．
逆に言えば，ARIMAモデルは指数平滑法よりも一般化されたモデルである．
ETSモデルがARIMAモデルでどの様に表現されるかは参考文献[11]が詳しい．
わざわざ特別解を求めても仕方が無いので，略．
－ARIMA（自己回帰和分移動平均）モデルによる予測（前回記事[1]の話）
-6500位のスコア
library(forecast)
x <- read.csv("data_6500.csv", header=F)
x <- x^T
xts <- ts(as.numeric(x), start=c(20, 16), frequency=24)
xts.stl <- stl(xts, s.window="periodic")
plot(xts.stl)
pp.test(xts)
library(tseries)
adf.test(xts)
library(forecast)
model <- auto.arima(xts, trace=T, stepwise=F, approximation=F)
summary(model)
tsdiag(model, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
plot(forecast(model, level=c(95,99.99),h=200), bty="n", pos=0)
axis(side=4, pos=20.625)
結果：
> xts <- ts(as.numeric(x), start=c(20, 16), frequency=24) > xts.stl <- stl(xts, s.window="periodic") > plot(xts.stl)
周期変動がみられるが，残差と比較してもスケールが小さく明確で無い．
> PP.test(xts) Phillips-Perron Unit Root Test data: xts Dickey-Fuller Z(alpha) = -18.6557, Truncation lag parameter = 3, p-value = 0.07546 alternative hypothesis: stationary > adf.test(xts) Augmented Dickey-Fuller Test data: xts Dickey-Fuller = -3.2457, Lag order = 4, p-value = 0.0871 alternative hypothesis: stationary
P値は小さいが棄却される程ではない．
> model <- auto.arima(xts, trace=T, stepwise=F, approximation=F) > summary(model) Series: xts ARIMA(0,1,2)(1,0,0)[24] with drift Coefficients: ma1 ma2 sar1 drift 1.1889 0.9069 0.3009 161.6462 s.e. 0.0587 0.0909 0.1611 22.4565 sigma^2 estimated as 2480: log likelihood=-391.59 AIC=793.18 AICc=794.07 BIC=804.7 Training set error measures: ME RMSE MAE MPE MAPE MASE Training set -4.568079 49.12618 32.16711 -Inf Inf 0.01188655 > > tsdiag(model, gof.lag=24)#残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット >
残差に有意な自己相関がみられる為，あまり良いモデルではない．
30.45833 37044.51 31884.78 42204.23 26802.29 47286.72
30.50000 37206.16 32028.38 42383.95 26928.10 47484.23
30.54167 37367.82 32172.05 42563.60 27054.04 47681.61
30.58333 37529.49 32315.77 42743.20 27180.10 47878.88
30.62500 37691.11 32459.52 42922.70 27306.24 48075.98
予測値：
37,691.11±5,231.59（95%CI）
37,691.11±10,384.87（99.99%CI）
－32,500位のスコア
x2 <- read.csv("data_32500.csv", header=F)
x2 <- x2^T
x2ts <- ts(as.numeric(x2), start=c(20, 16), frequency=24)
x2ts.stl <- stl(x2ts, s.window="periodic")
plot(x2ts.stl)
PP.test(x2ts)
library(tseries)
adf.test(x2ts)
library(forecast)
model2 <- auto.arima(x2ts, trace=T, stepwise=F, approximation=F)
summary(model2)
tsdiag(model2, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
plot(forecast(model2, level=c(95,99.99),h=200), bty="n", pos=0)
axis(side=4, pos=20.625)
結果：
x2ts <- ts(as.numeric(x2), start=c(20, 16), frequency=24) x2ts.stl <- stl(x2ts, s.window="periodic") plot(x2ts.stl)
> PP.test(x2ts) Phillips-Perron Unit Root Test data: x2ts Dickey-Fuller Z(alpha) = -13.2646, Truncation lag parameter = 3, p-value = 0.3309 alternative hypothesis: stationary > adf.test(x2ts) Augmented Dickey-Fuller Test data: x2ts Dickey-Fuller = -2.8587, Lag order = 4, p-value = 0.2252 alternative hypothesis: stationary
検定結果から本データは単位根過程に従うと考えて良い．
> model2 <- auto.arima(x2ts, trace=T, stepwise=F, approximation=F) > summary(model2) Series: x2ts ARIMA(4,1,0) with drift Coefficients: ar1 ar2 ar3 ar4 drift 0.7067 -0.4395 0.9482 -0.6028 67.9292 s.e. 0.1251 0.1217 0.0781 0.1105 7.2805 sigma^2 estimated as 560.7: log likelihood=-336.74 AIC=685.48 AICc=686.73 BIC=699.31 Training set error measures: ME RMSE MAE MPE MAPE MASE Training set -0.6401462 23.36032 17.18989 -Inf Inf 0.01444856
> tsdiag(model2, gof.lag=24)
残差に有意な自己相関はみられず，妥当なモデルであると考えられる．
> plot(forecast(model2, level=c(95,99.99),h=200), bty="n", pos=0) > axis(side=4, pos=20.625)
30.41667 15808.003 14290.543 17325.463 12795.796 18820.211
30.45833 15876.107 14353.991 17398.222 12854.658 18897.555
30.50000 15944.156 14417.325 17470.986 12913.348 18974.964
30.54167 16011.823 14480.246 17543.401 12971.593 19052.054
30.58333 16079.840 14543.637 17616.042 13030.428 19129.252
30.62500 16147.954 14607.102 17688.806 13089.312 19206.596
予測値：
16,148±1,541（95%CI）
16,148±3,059（99.99%CI）
－ARFIMA（自己回帰実数和分移動平均）モデルによる予測（library(fracdiff)，fdGPH()やfracdiff()）
－GARCH（一般化自己回帰条件付分散不均一）モデルによる予測（library(fSeries)，garchfit()関数）
－SETAR（自己励起型閾値）モデルによる予測（library(tsDyn)，setar()関数）
は本データではやる意味が無いので略．ARIMAで十分当て嵌まりが良い．
重要なのはこんだけ当て嵌まりが良くて，予測もうまく出来ていそうなのに，
インパクトがあると（それまでと全く異なる挙動が含まれると），一気に誤差が生じるという事．
（要はEX解禁と共に一気にスコアが跳ね上がるのは分かり切っているけど，
それと共にこれらのモデルが説明力を失う事も分かり切っているという事……）
その部分をどの様に考慮するかが大事な話なんだけど，また，後日続きを書く．
●2014/6/20（金） 16:00～2014/6/28（土） 18:35のデータから予測
－線形回帰モデルによる予測（lm()関数）
－一般化線形混合モデルによる予測（library(nlme)，lme()関数）
－ARIMA（自己回帰和分移動平均）モデルによる予測（前回記事[1]の話）
－ARFIMA（自己回帰実数和分移動平均）モデルによる予測（library(fracdiff)，fdGPH()やfracdiff()）
－GARCH（一般化自己回帰条件付分散不均一）モデルによる予測（library(fSeries)，garchfit()関数）
－SETAR（自己励起型閾値）モデルによる予測（library(tsDyn)，setar()関数）
を考える．
－6500位のボーダー予測
・線形回帰モデルによる予測
xy <- read.csv("data_6500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
xy_lmmodel <- lm(y~x-1, xy)#単回帰（線形回帰）モデル：y=ax（切片は0）
summary(xy_lmmodel)
par(mfrow=c(2,2))
plot(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel, se.fit=T, interval="confidence")
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_lmpred_se <- 2*(sefit$coefficients[2])
new_x <- seq(0, 240,1)
new_y <- new_x*xy_lmmodel_coefficients#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_lmpred_se
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_lmmodel) Call: lm(formula = y ~ x - 1, data = xy) Residuals: Min 1Q Median 3Q Max -2452.3 -531.8 701.7 1098.2 1852.4 Coefficients: Estimate Std. Error t value Pr(>|t|) x 131.0628 0.7667 170.9 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1201 on 194 degrees of freedom Multiple R-squared: 0.9934, Adjusted R-squared: 0.9934 F-statistic: 2.922e+04 on 1 and 194 DF, p-value: < 2.2e-16
> par(mfrow=c(2,2)) > plot(xy_lmmodel)
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
[236,] NA 30799.7596 30675.91230 30923.60694
[237,] NA 30930.8224 30806.49632 31055.14853
[238,] NA 31061.8852 30937.08034 31186.69013
[239,] NA 31192.9480 31067.66436 31318.23172
[240,] NA 31324.0108 31198.24838 31449.77331
[241,] NA 31455.0737 31328.83240 31581.31491
予測値：31455.0737±123.24121（95%CI）
これはあくまで母平均推定値．これにモデル残差を加えて勘案すると，
今後，残差は増大する事が予想される為，およそ20毎に残差2,000とすると
予測ボーダー（目安）：約3,8000 安全圏（倍値）：43,000
－一般化線形混合モデルによる予測（library(nlme)，lme()関数）
Rコマンド
library(nlme)
xy <- read.csv("data_6500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
atemp <- rep(1:24, 10)
a <- atemp[1:length(xy[,1])]
xy <- cbind(xy, a)
gxy <-groupedData(y ~ x | a, data =xy)
xy_lmmodel <- lme(y~x, gxy)
summary(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel)
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
xy_coeff <- xy_lmmodel_coefficients$fixed #[1]が切片，[2]がxの係数
xy_random <- xy_lmmodel_coefficients$random #$a[a,1]が切片，$a[a,2]が係数
new_x <- seq(0, 240,1)
new_a <- rep(1:24,11)
new_a <- new_a[1:length(new_x)]
new_y <- xy_coeff[1] + new_x*xy_coeff[2] + xy_random$a[new_a,1] + new_x*xy_random$a[new_a, 2]#将来の値を予測
new_CI <- cbind(new_y-2*(sd(DATA[,2])/sqrt(length(xy[,2]))), new_y+2*(sd(DATA[,2])/sqrt(length(xy[,2]))))#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
plot(xy_lmmodel)
> summary(xy_lmmodel) Linear mixed-effects model fit by REML Data: gxy AIC BIC logLik 3273.434 3293.01 -1630.717 Random effects: Formula: ~x | a Structure: General positive-definite StdDev Corr (Intercept) 7.550795e-02 (Intr) x 6.927583e-04 -0.864 Residual 1.077252e+03 Fixed effects: y ~ x Value Std.Error DF t-value p-value (Intercept) 1065.4042 153.69566 170 6.93191 0 x 122.8463 1.37044 170 89.63974 0 Correlation: (Intr) x -0.865 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.2586355 -0.6787230 0.3021622 0.7115715 1.6388087 Number of Observations: 195 Number of Groups: 24
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
23 NA 29934.290 29294.9865 30573.593
24 NA 30057.136 29417.8329 30696.440
1 NA 30179.983 29540.6792 30819.286
2 NA 30302.829 29663.5255 30942.132
3 NA 30425.675 29786.3718 31064.979
4 NA 30548.521 29909.2181 31187.825
予測値：30,548±639.304
－ARIMA（自己回帰和分移動平均）モデルによる予測（前回記事[1]の話）
Rコマンド
library(forecast)
x <- read.csv("data_6500.csv", header=F)
x <- x^T
xts <- ts(as.numeric(x), start=c(20, 16), frequency=24)
xts.stl <- stl(xts, s.window="periodic")
plot(xts.stl)
PP.test(xts)
adf.test(xts)
model <- auto.arima(xts, trace=T, stepwise=F, approximation=F)
summary(model)
tsdiag(model, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
plot(forecast(model, level=c(95,99.99),h=200), bty="n", pos=0)
axis(side=4, pos=20.625)
結果：
> xts <- ts(as.numeric(x), start=c(20, 16), frequency=24) > xts.stl <- stl(xts, s.window="periodic") > plot(xts.stl)
心電図の様な周期変動がみられる．
前半と後半で綺麗な線形トレンドがみられる（26日でデータを区切って別に考えた方が，当て嵌まりの良いモデルを立てられそう）．
> PP.test(xts) Phillips-Perron Unit Root Test data: xts Dickey-Fuller = -0.6554, Truncation lag parameter = 4, p-value = 0.973 > adf.test(xts) Augmented Dickey-Fuller Test data: xts Dickey-Fuller = -1.0642, Lag order = 5, p-value = 0.9258 alternative hypothesis: stationary
> model <- auto.arima(xts, trace=T, stepwise=F, approximation=F) > summary(model) Series: xts ARIMA(1,1,3)(1,0,0)[24] with drift Coefficients: ar1 ma1 ma2 ma3 sar1 drift 0.5843 0.6587 0.6065 0.3590 0.2307 153.0758 s.e. 0.1189 0.1197 0.1320 0.1233 0.0854 24.0923 sigma^2 estimated as 1793: log likelihood=-998.12 AIC=2010.23 AICc=2010.84 BIC=2033.11 Training set error measures: ME RMSE MAE MPE MAPE MASE Training set -2.39467 42.12178 28.58586 -Inf Inf 0.005063638 > tsdiag(model, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
> plot(forecast(model, level=c(95,99.99),h=50), bty="n", pos=0) > axis(side=4, pos=20.625)
30.41667 33132.96 29629.31 36636.61 26178.11 40087.81
30.45833 33286.30 29723.84 36848.75 26214.71 40357.89
30.50000 33440.11 29819.80 37060.43 26253.68 40626.55
30.54167 33594.09 29916.83 37271.35 26294.61 40893.57
30.58333 33749.56 30016.22 37482.89 26338.76 41160.35
30.62500 33902.79 30114.20 37691.37 26382.33 41423.25
予測値：
33,902.79±3788.58（95%CI）
33,902.79±7520.21（99.99%CI）
予測ボーダー（目標）：42,000（99.99%CI上限程度を目安）
予測ボーダー（安全圏）：45268.53（シックスシグマ）
－ARFIMA（自己回帰実数和分移動平均）モデルによる予測（library(fracdiff)，fdGPH()やfracdiff()）
書きかけ
－32,500位のボーダー予測
・線形回帰モデルによる予測
まずは普通に線形回帰してみる．
xy <- read.csv("data_32500_2.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
xy_lmmodel <- lm(y~x-1, xy)#単回帰（線形回帰）モデル：y=ax（切片は0）
summary(xy_lmmodel)
par(mfrow=c(2,2))
plot(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel, se.fit=T, interval="confidence")
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_lmpred_se <- 2*(sefit$coefficients[2])
new_x <- seq(0, 240,1)
new_y <- new_x*xy_lmmodel_coefficients#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_lmpred_se
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
結果：
> summary(xy_lmmodel) Call: lm(formula = y ~ x - 1, data = xy) Residuals: Min 1Q Median 3Q Max -1548.5 -917.6 -297.2 182.5 2265.0 Coefficients: Estimate Std. Error t value Pr(>|t|) x 76.1463 0.5198 146.5 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 929.6 on 212 degrees of freedom Multiple R-squared: 0.9902, Adjusted R-squared: 0.9902 F-statistic: 2.146e+04 on 1 and 212 DF, p-value: < 2.2e-16
> par(mfrow=c(2,2)) > plot(xy_lmmodel)
> matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
見るからに当て嵌まりは宜しくないし，残差も増大している事が明らかなので，
改めてデータを区切って考えてみる．残差のグラフをみると，
スコア10,000（t=150）の所からトレンドがみられる．
念の為に確認してみると，
#stlは時系列データ(Time Series)を与えないといけないので改めてデータ作成
> x2 <- read.csv("data_32500.csv", header=F)#xy[,2]（スコア値）に同じ
> x2 <- x2^T
> x2ts <- ts(as.numeric(x2), start=c(20, 16), frequency=24)
> x2ts.stl <- stl(x2ts, s.window="periodic")
> plot(x2ts.stl)
20（金）～26日（木）まで綺麗なトレンドがみられ，26日に増大，その後は少し緩やかになるも，
それまでよりは強いupward biasをもったトレンドがみられる．
これらを勘案すると，27日以降のデータで線形回帰を行ない，
上振れ値として26日の増分を加味すれば，
割と妥当な予測値が得られるのではないだろうか．
簡単な手法（単純な線形回帰）と簡単な仮定を置くだけで，
それなりに当て嵌まりの良い目標値が得られれば嬉しいよね．
26日の増分：3,406（ちなみに6,500位の場合は5,751）
27日以降のトレンド（線形回帰）は
xy <- read.csv("data_32500_3.csv")#説明変数（独立変数）と目的変数（従属変数）の構造のデータセットに
a <- 11741 #06/27（金）00:35の値（開始値）
xy_lmmodel <- lm(y-a~x-1, xy)#単回帰（線形回帰）モデル：y=ax（切片は0）を求める
#データ全体から開始時の値を引いて切片0の回帰式を求めてから後で切片としてaを足す
summary(xy_lmmodel)
xy_lmpred <- predict(xy_lmmodel, se.fit=T, interval="confidence")
xy_lmmodel_coefficients <- xy_lmmodel$coefficients#y=axの係数aを取り出す
retu <- 1:length(xy[,1])
sedata <- xy_glmpred$se.fit
DF <- data.frame(N=retu, SE=sedata)
sefit <- lm(SE~N, DF)
xy_lmpred_se <- 2*(sefit$coefficients[2])
new_x <- seq(0, 240,1)
new_y <- a+new_x*xy_lmmodel_coefficients#将来の値を予測
new_yse <- sefit$coefficients[1]+new_x*xy_lmpred_se
new_CI <- cbind(new_y-new_yse, new_y+new_yse)#95CIのデータ系列
DATA <- cbind(xy[,2], new_y, new_CI)#プロット用のデータセットを作成
DATA[, 1] <- replace(DATA[,1], (length(xy[,1])+1):241 , NA)#ゴミはNAに
matplot(new_x, DATA, lty=c(1,2,3,3), type="l", bty="n", xlab="", ylab="", xlim=c(0, 250), ylim=c(0, 40000), pos=0)
> par(mfrow=c(2,2)) > plot(xy_lmmodel)
[84,] NA 20289.66 20238.59 20340.73
[85,] NA 20392.66 20341.11 20444.21
[86,] NA 20495.65 20443.63 20547.68
[87,] NA 20598.65 20546.14 20651.16
[88,] NA 20701.65 20648.66 20754.63
予測値：20,701.65±52.98
予想ボーダー（目標値）：21,901（割と妥当な値だと思うが果たして）
予想ボーダー（安全圏）：24,101（今のトレンドに26日並みのインパクトが加算されたとして）
－ARIMA（自己回帰和分移動平均）モデルによる予測（前回記事[1]の話）
Rコマンド
x2 <- read.csv("data_32500.csv", header=F)
x2 <- x2^T
x2ts <- ts(as.numeric(x2), start=c(20, 16), frequency=24)
x2ts.stl <- stl(x2ts, s.window="periodic")
plot(x2ts.stl)
PP.test(x2ts)
library(tseries)
adf.test(x2ts)
library(forecast)
model2 <- auto.arima(x2ts, trace=T, stepwise=F, approximation=F)
summary(model2)
tsdiag(model2, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
plot(forecast(model2, level=c(95,99.99),h=200), bty="n", pos=0)
axis(side=4, pos=20.625)
結果：
> x2ts <- ts(as.numeric(x2), start=c(20, 16), frequency=24) > x2ts.stl <- stl(x2ts, s.window="periodic") > plot(x2ts.stl)
> PP.test(x2ts) Phillips-Perron Unit Root Test data: x2ts Dickey-Fuller = 0.7367, Truncation lag parameter = 4, p-value = 0.99 > adf.test(x2ts) Augmented Dickey-Fuller Test data: x2ts Dickey-Fuller = 0.0295, Lag order = 5, p-value = 0.99 alternative hypothesis: stationary 警告メッセージ: In adf.test(x2ts) : p-value greater than printed p-value #「p値クッソ大きいぞほげえええ」って怒られているけど #帰無仮説（単位根過程である）が棄却されない事を確認しているので気にしない
> model2 <- auto.arima(x2ts, trace=T, stepwise=F, approximation=F) > summary(model2) Series: x2ts ARIMA(4,2,0)(1,0,0)[24] Coefficients: ar1 ar2 ar3 ar4 sar1 -0.1157 -0.3954 0.1511 -0.2240 0.3574 s.e. 0.0790 0.0752 0.0765 0.0789 0.0785 sigma^2 estimated as 881.2: log likelihood=-1016.78 AIC=2045.55 AICc=2045.96 BIC=2065.66 Training set error measures: ME RMSE MAE MPE MAPE MASE Training set -0.2784047 29.54532 18.72698 -0.3368363 0.7845073 0.004288119 > tsdiag(model2, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
> plot(forecast(model2, level=c(95,99.99),h=50), bty="n", pos=0) > axis(side=4, pos=20.625)
30.41667 22369.70 19821.99 24917.41 17312.41 27426.98
30.45833 22558.02 19851.74 25264.30 17185.97 27930.07
30.50000 22743.79 19868.80 25618.77 17036.85 28450.72
30.54167 22928.50 19876.11 25980.88 16869.41 28987.59
30.58333 23114.51 19879.16 26349.87 16692.23 29536.80
30.62500 23299.29 19873.98 26724.59 16499.94 30098.63
予測値：
23299.29±3425.3（95%CI）
23299.29±6799.34（99.99%CI）
> par(mfrow=c(3,1)) > plot(x2ts) > plot(diff(x2ts)) > plot(diff(diff(x2ts))) > abline(h=0)
26日のインパクトがでか過ぎてそれに引っ張られているのがよく分かる．
（後述するがこの部分を別の挙動と捉えてうまく考えようとするのがGARCHやSETAR）
それ以外の部分をみると，1回階差より2回階差の方が規則的な周期変動のグラフになっている．
この様に（インパクトを除けば），一般的なデータは階差を取れば取る程
当て嵌まりが良くなるので，高階差のモデリングをしてしまうリスクがある．
それを改善したモデルがARFIMA（自己回帰実数和分移動平均）モデル．
今のデータの場合，元データは単位根過程に従うが，1階階差をとれば
統計的には定常過程に従うと考えられるので，1階階差でモデリングし直してみる．
（どちらにせよ，ARIMAではインパクトによる誤差はどうしようもないが）
> model2 <- auto.arima(x2ts, d=1,trace=T, stepwise=F, approximation=F) > summary(model2) Series: x2ts ARIMA(4,1,0)(1,0,0)[24] with drift Coefficients: ar1 ar2 ar3 ar4 sar1 drift 0.7988 -0.2489 0.4823 -0.2443 0.3554 93.8348 s.e. 0.0782 0.0964 0.0954 0.0774 0.0769 13.7909 sigma^2 estimated as 854.1: log likelihood=-1013.32 AIC=2040.65 AICc=2041.2 BIC=2064.14 Training set error measures: ME RMSE MAE MPE MAPE MASE Training set -0.7169863 29.08824 18.62839 -Inf Inf 0.004265543 > tsdiag(model2, gof.lag=24) #残差分析，上から残差，残差自己相関，Ljung-Box検定のp値のプロット
> plot(forecast(model2, level=c(95,99.99),h=200), bty="n", pos=0) > axis(side=4, pos=20.625)
30.45833 21309.57 20150.08 22469.06 19007.94 23611.19
30.50000 21429.19 20233.92 22624.45 19056.55 23801.83
30.54167 21544.95 20311.00 22778.91 19095.52 23994.39
30.58333 21659.20 20385.78 22932.61 19131.43 24186.96
30.62500 21769.09 20454.44 23083.73 19159.47 24378.70
予測値：
21769.09±1314.64（95%CI）
21769.09±2609.61（99.99%CI）
予想ボーダー（目標値）：23,083～23,300（95%CI上限～2階階差の結果程度が目安か）
予想ボーダー（安全圏）：25,713（シックスシグマ）
参考文献
[1] ソシャゲ（スクフェス）のボーダーを統計的に予測する方法-非定常時系列データ分析-｜粉末@それは風のように (日記)
[2] 沖本竜義：経済・ファイナンスデータの計量時系列分析．朝倉書店，2010/02．
[3] 金 明哲：Rによるデータサイエンス – データ解析の基礎から最新手法まで．森北出版，2007/10．
[4] 辻谷 将明，竹澤 邦夫，金 明哲：マシンラーニング Rで学ぶデータサイエンス 6．共立出版，2009/06．
[5] 統計解析フリーソフト R の備忘録頁 ver.3.1｜統計解析編｜R-tips
[6] Rと非線形回帰分析｜[連載] フリーソフトによるデータ解析・マイニング 第16回
[7] 平滑化スプラインと加法モデル｜海と魚と統計解析
[8] Rで計量時系列分析～CRANパッケージ総ざらい～｜slideshare.net
[9] Rで計量時系列分析：状態変化を伴うモデル（閾値モデル、平滑推移モデル、マルコフ転換モデル）｜銀座で働くデータサイエンティストのブログ
[10] R for MATLAB users
[11] 8.10 ARIMA vs ETS｜OTexts
[12] 生態学のデータ解析 – Gamma 分布の GLM｜生態学のデータ解析
粉末@それは風のように (日記) 