そう言えば、ちょっと前のデータ分析業界5年間振り返り記事で「人工知能ブームに引っ張られてデータサイエンティストブームも再燃しつつある」みたいなことを書いたわけですが、本当にそうなんだっけ?というところをこれまでに検証したことはなかったなぁと思い出したのでした。
ということで、手っ取り早く計量時系列分析の手法を使ってこの2つのブームの関係性を解明してみようと思います。ちなみにデータソースは以下の2つのGoogleトレンド検索結果です。それぞれ'multiTimeline_ai.csv', 'multiTimeline_ds.csv'みたいなファイル名にして保存してあります。
やることは基本的には完全に決まっていますので、お定まりのルーチンをただこなすだけです。では、早速やってみましょう。
データを読み込む
上記のようにCSVファイルで保存してあるはずなので、ベタッと読み込んでついでにプロットしてみます。
> d1 <- read.csv('multiTimeline_ai.csv', skip = 2) > d2 <- read.csv('multiTimeline_ds.csv', skip = 2) > d <- data.frame(ai=d1[,2], ds=d2[,2]) > head(d) ai ds 1 3 2 2 9 9 3 5 13 4 7 22 5 7 28 6 3 28 > matplot(d, type='l', col=c(2,4), lwd=3, lty=1, ylab='') > legend('bottomright', col=c(2,4), lwd=3, lty=1, legend=c('人工知能', 'データサイエンティスト'), ncol=2)
これでこの記事のトップと同じプロットが得られたはずです。見た感じでは、お互いに微妙な均衡関係にあるような気もしますが。。。という話題は後できっちりやるとして、ここでは一旦全体の傾向を確認するに留めておきます。原系列はweeklyで尚且つ特に周期性のないデータのはずですが、念のためFFTでパワースペクトルを見ておくことにしましょう。
> par(mfrow=c(2,1)) > plot(abs(fft(d$ai)), type='l', lwd=2, xlab='', ylab='', main='AI') > plot(abs(fft(d$ds)), type='l', lwd=2, xlab='', ylab='', main='DS')
ぱっと見では特徴的な周波数は(低周波領域を除いて)なさそうです。よって、特に季節調整は考慮せずに以下の分析に移ることとします。
単位根過程まわりのチェックをする
見た感じでは長期トレンドがあるようにも見える2つの時系列なので、一応ADF検定してみます。
> library(tseries) > adf.test(d$ai) Augmented Dickey-Fuller Test data: d$ai Dickey-Fuller = -3.4571, Lag order = 6, p-value = 0.04736 alternative hypothesis: stationary > adf.test(d$ds) Augmented Dickey-Fuller Test data: d$ds Dickey-Fuller = -3.8939, Lag order = 6, p-value = 0.01479 alternative hypothesis: stationary
p値至上主義を信じているわけではないのですが。。。ギリギリどちらも「定常過程」という判定になりました。が、これだけだと心許ないのでKPSS検定、PP検定も試してみます。
> pp.test(d$ai) Phillips-Perron Unit Root Test data: d$ai Dickey-Fuller Z(alpha) = -85.993, Truncation lag parameter = 5, p-value = 0.01 alternative hypothesis: stationary 警告メッセージ: pp.test(d$ai) で: p-value smaller than printed p-value > pp.test(d$ds) Phillips-Perron Unit Root Test data: d$ds Dickey-Fuller Z(alpha) = -152.43, Truncation lag parameter = 5, p-value = 0.01 alternative hypothesis: stationary 警告メッセージ: pp.test(d$ds) で: p-value smaller than printed p-value > library(urca) > summary(ur.kpss(d$ai)) ####################### # KPSS Unit Root Test # ####################### Test is of type: mu with 5 lags. Value of test-statistic is: 3.7547 Critical value for a significance level of: 10pct 5pct 2.5pct 1pct critical values 0.347 0.463 0.574 0.739 > summary(ur.kpss(d$ds)) ####################### # KPSS Unit Root Test # ####################### Test is of type: mu with 5 lags. Value of test-statistic is: 0.5191 Critical value for a significance level of: 10pct 5pct 2.5pct 1pct critical values 0.347 0.463 0.574 0.739
うーむ。。。KPSS検定の場合に限って「人工知能」「データサイエンティスト」とも単位根過程という結果になりました。この後の因果性を見る際は、念のためケースバイケースながら差分系列に直すか、共和分関係をチェックした方が良さそうです。
VARモデルを推定し、共和分をチェックする
とるものもとりあえず、2変数VARモデルを原系列・差分系列それぞれについて推定します。
> library(vars) > VARselect(d, lag.max=52, type='both') $selection AIC(n) HQ(n) SC(n) FPE(n) 3 1 1 3 $criteria 1 2 3 4 5 6 7 AIC(n) 8.783978 8.781656 8.769103 8.788558 8.779134 8.806056 8.808841 HQ(n) 8.835883 8.859513 8.872913 8.918320 8.934848 8.987724 9.016461 SC(n) 8.912345 8.974206 9.025837 9.109475 9.164234 9.255340 9.322308 FPE(n) 6528.857819 6513.861234 6432.890664 6559.739234 6498.912587 6677.265659 6697.233792 8 9 10 11 12 13 14 AIC(n) 8.809953 8.836354 8.861409 8.889940 8.900645 8.915896 8.913674 HQ(n) 9.043526 9.095879 9.146887 9.201370 9.238027 9.279231 9.302961 SC(n) 9.387604 9.478189 9.567427 9.660141 9.735029 9.814464 9.876425 FPE(n) 6706.424453 6888.071494 7065.632462 7273.586763 7356.015645 7473.990811 7463.089330 15 16 17 18 19 20 21 AIC(n) 8.899206 8.928812 8.950994 8.979221 9.005974 9.014756 9.019717 HQ(n) 9.314446 9.370004 9.418139 9.472319 9.525024 9.559758 9.590672 SC(n) 9.926141 10.019930 10.106296 10.198706 10.289643 10.362608 10.431752 FPE(n) 7362.325509 7591.091872 7770.065443 8002.572592 8231.084389 8316.612440 8372.339266 22 23 24 25 26 27 28 AIC(n) 9.023659 9.009730 9.007327 9.024690 9.043345 9.049487 9.082855 HQ(n) 9.620566 9.632589 9.656139 9.699455 9.744062 9.776157 9.835477 SC(n) 10.499878 10.550132 10.611913 10.693459 10.776297 10.846623 10.944174 FPE(n) 8421.301543 8322.019869 8320.835492 8487.428595 8670.374564 8749.056759 9074.215962 29 30 31 32 33 34 35 AIC(n) 9.112045 9.136369 9.160659 9.182594 9.213515 9.245259 9.241117 HQ(n) 9.890619 9.940896 9.991138 10.039026 10.095900 10.153595 10.175406 SC(n) 11.037548 11.126056 11.214529 11.300647 11.395752 11.491678 11.551720 FPE(n) 9374.496745 9640.154883 9915.626476 10177.811794 10544.292149 10936.265305 10946.452431 36 37 38 39 40 41 42 AIC(n) 9.236669 9.267494 9.242443 9.264997 9.272241 9.283538 9.315826 HQ(n) 10.196911 10.253689 10.254590 10.303096 10.336292 10.373542 10.431783 SC(n) 11.611456 11.706464 11.745597 11.832334 11.903761 11.979242 12.075713 FPE(n) 10956.924010 11365.075980 11151.822492 11480.391606 11643.643709 11862.033571 12346.124028 43 44 45 46 47 48 49 AIC(n) 9.320626 9.287883 9.281527 9.288248 9.30243 9.311838 9.330188 HQ(n) 10.462535 10.455744 10.475341 10.508015 10.54815 10.583510 10.627812 SC(n) 12.144697 12.176137 12.233964 12.304869 12.38323 12.456826 12.539359 FPE(n) 12507.127547 12209.029734 12242.644852 12444.235005 12750.60147 13009.462149 13400.510196 50 51 52 AIC(n) 9.333244 9.302537 9.296769 HQ(n) 10.656820 10.652066 10.672251 SC(n) 12.606598 12.640075 12.698491 FPE(n) 13601.992099 13356.545984 13455.567081 > d.var <- VAR(d, p=3, type='both') > dd <- data.frame(ai=diff(d$ai), ds=diff(d$ds)) > VARselect(dd, lag.max=52, type='const') $selection AIC(n) HQ(n) SC(n) FPE(n) 7 4 2 7 $criteria 1 2 3 4 5 6 7 AIC(n) 9.069478 8.967836 8.941902 8.889113 8.884470 8.847570 8.832094 HQ(n) 9.108543 9.032944 9.033053 9.006307 9.027706 9.016849 9.027416 SC(n) 9.166079 9.128837 9.167304 9.178915 9.238672 9.266172 9.315097 FPE(n) 8686.126028 7846.754545 7646.122412 7253.382060 7220.431941 6959.763286 6854.103755 8 9 10 11 12 13 14 AIC(n) 8.856046 8.891046 8.909209 8.918474 8.932824 8.919231 8.902155 HQ(n) 9.077411 9.138454 9.182660 9.217968 9.258361 9.270811 9.279778 SC(n) 9.403449 9.502849 9.585413 9.659078 9.737829 9.788636 9.835960 FPE(n) 7021.889685 7274.147208 7410.181063 7482.457174 7594.608001 7496.718046 7375.100777 15 16 17 18 19 20 21 AIC(n) 8.932496 8.960475 8.989160 9.014198 9.023604 9.034094 9.033683 HQ(n) 9.336162 9.390183 9.444911 9.495993 9.531441 9.567974 9.593607 SC(n) 9.930702 10.023081 10.116166 10.205605 10.279411 10.354301 10.418291 FPE(n) 7608.611519 7831.903334 8068.458355 8283.011847 8372.568183 8473.554413 8484.145661 22 23 24 25 26 27 28 AIC(n) 9.024043 9.030691 9.049486 9.063996 9.074963 9.104839 9.137462 HQ(n) 9.610009 9.642700 9.687538 9.728092 9.765101 9.821020 9.879686 SC(n) 10.473051 10.544100 10.627295 10.706206 10.781572 10.875849 10.972872 FPE(n) 8418.131509 8491.317332 8671.449597 8819.287234 8939.779122 9236.935766 9572.410413 29 30 31 32 33 34 35 AIC(n) 9.159901 9.189402 9.206854 9.236028 9.270028 9.256824 9.261208 HQ(n) 9.928168 9.983712 10.027207 10.082424 10.142466 10.155306 10.185732 SC(n) 11.059712 11.153613 11.235466 11.329040 11.427440 11.478637 11.547421 FPE(n) 9821.925183 10151.925508 10370.082480 10720.781590 11140.173603 11045.577519 11149.605599 36 37 38 39 40 41 42 AIC(n) 9.284935 9.274826 9.294455 9.313068 9.321081 9.355579 9.371389 HQ(n) 10.235502 10.251437 10.297108 10.341764 10.375820 10.436361 10.478214 SC(n) 11.635548 11.689840 11.773869 11.856882 11.929296 12.028194 12.108404 FPE(n) 11478.216457 11427.307834 11724.212974 12021.126252 12200.358246 12719.768775 13021.305246 43 44 45 46 47 48 49 AIC(n) 9.357677 9.37988 9.383846 9.390585 9.412743 9.45004 9.454264 HQ(n) 10.490545 10.53879 10.568800 10.601582 10.649783 10.71312 10.743389 SC(n) 12.159093 12.24570 12.314063 12.385202 12.471760 12.57346 12.642082 FPE(n) 12947.963597 13352.01721 13526.403399 13748.037148 14197.913520 14894.37571 15125.357388 50 51 52 AIC(n) 9.430069 9.376903 9.396024 HQ(n) 10.745238 10.718115 10.763279 SC(n) 12.682287 12.693522 12.777043 FPE(n) 14938.538524 14341.711746 14810.613641 > dd.var <- VAR(dd, p=7, type='const')
先にも書いたように、原系列は単位根過程同士である可能性があります。ということで、Johansenの方法で共和分検定を行います。
> d.vecm <- ca.jo(d, ecdet='trend', type='eigen', K=3, spec='transitory') > summary(d.vecm) ###################### # Johansen-Procedure # ###################### Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration Eigenvalues (lambda): [1] 1.077463e-01 1.068366e-01 1.387779e-17 Values of teststatistic and critical values of test: test 10pct 5pct 1pct r <= 1 | 29.04 10.49 12.25 16.26 r = 0 | 29.30 16.85 18.96 23.65 Eigenvectors, normalised to first column: (These are the cointegration relations) ai.l1 ds.l1 trend.l1 ai.l1 1.0000000 1.0000000 1.000000 ds.l1 13.1803243 -0.1591407 -3.051548 trend.l1 -0.3832159 -0.2174744 25.629559 Weights W: (This is the loading matrix) ai.l1 ds.l1 trend.l1 ai.d -0.0008666047 -0.30194326 -3.340989e-19 ds.d -0.0276692419 -0.02267098 2.786357e-20
ガッチガチに共和分が入っていますorz ということで、ベクトル誤差修正モデルに変換します。またこの結果とGranger表現定理を踏まえると、差分系列のVARモデルは存在し得ない(つまり推定しても正しいモデルではない=dd.varはモデルとして不適切)ということになります。
> d.vec2var <- vec2var(d.vecm, r=1)
これで準備が整いました。ちなみに共和分関係にあるということは、沖本本の記述に従えば「人工知能ブームが過熱するとデータサイエンティストブームが萎む」という均衡関係にある可能性が示唆されるということでもあります。
未来値の予測と、因果性の判定を行う
未来値の予測はベクトル誤差修正モデルにpredictメソッドをかければおしまいです。
> plot(predict(d.vec2var, n.ahead=26, ci=0.95), lwd=2)
「人工知能」はわずかながら上昇トレンドを維持する一方で、「データサイエンティスト」は今後伸び悩むという予測結果になっています。因果性の判定はベクトル誤差修正モデルだとGranger因果を使えないので、インパルス応答関数推定でやってみましょう。
> plot(irf(d.vec2var, impulse='ds', response = 'ai', n.ahead=26, ci=0.95), lwd=2) > plot(irf(d.vec2var, impulse='ai', response = 'ds', n.ahead=26, ci=0.95), lwd=2)
どうやら因果性はなさそうだという結果になりました。