注意: この記事のRパッケージ解説は公式ドキュメントではありません. パッケージは常に改良されているため, なるべく最新の公式ドキュメントを参考にしてください.
概要
いまのところパッケージ側で修正される予定はなしだが, 放置すると忘れそうなのでとりあえず書いておく.
R-wakalang で, パネルデータを扱う plm パッケージが提供している cipstest 関数を使うと発生するエラーの解決法を知りたいという投稿があった. ソースコードを確認したところ, 投稿者が使い方を間違えているのではなく実装に問題があるとしか思えなかった. しかし適切な実装に直すには時間がかかる (し私が積極的に直しに行く予定はない) ため, 取り急ぎ周知してもらうだけにする.
なお, この問題を確認した plm パッケージのバージョンは 2.2-5 (R本体は 3.6.5 または 4.0.3) である.
問題
問題になった plm::cipstest 関数はヘルプドキュメントによれば Pesaran (2007) に基づいた cross-sectionally augmented IPS test for unit roots in panel models” を実装しているという. いわゆるパネルデータに対する単位根検定である.
しかし投稿者が使用したところ次のようなエラーが発生したという.
Error in cvals[nintl:ninth, tintl, i] : subscript out of bounds
日本語ロケールならば以下のようになるだろう.
cvals[nintl:ninth, tintl, i] でエラー: 添え字が許される範囲外です
このエラーの再現例も用意した.
require(plm) packageVersion("plm") ## 2.2.5 d <- expand.grid(i = 1:500, t = 1:20) set.seed(42) d$x1 <- rnorm(NROW(d)) d <- pdata.frame(d, index = c("i", "t")) pdim(d) ## Balanced Panel: n = 500, T = 20, N = 10000 cipstest(d$x1, type = "trend")
原因の詳細
以降は基本的に既に R-wakalang で回答したのと同じ話.
変数名からして内部処理の参照エラーっぽいのでソースコードを確認してみたところ, この関数の実装は, 「計算済みの数表から臨界値 (critical value) とp値を参照する」というものだった. そしてその表は, 提案元の Pesaran (2007) の論文に掲載されている数表の値そのままであるようだ.
この表はパネルデータの時系列の長さ () と個体数 (
) と有意水準の組み合わせごとに臨界値を掲載しているが,
それぞれ
までであり, 投稿者の使用したデータはもっと大きいものだったため, 参照エラーを起こした (200未満であれば表にない
の組み合わせに対しては線形補間した値を返している). 冒頭のエラーは
が大きい場合だったが,
が大きい場合はおそらく以下のようなエラーが出る.
cvals[which(rnam == n), tinth, i] でエラー: 添え字が許される範囲外です
つまりver. 2.2-5 現在, plm::cipstest 関数は または
が 200より大きなパネルデータでは使えないし, 返す値は厳密ではない.
解説
投稿者に上記の結果を回答したところ, 無視できる誤差だと判断して数表の中の最大の数値で代用することに妥協したようだ. 事実, 臨界値は単調減少するので, もしデータと目的に対して無視して問題ない誤差なら (例: 実際より小さいサイズでもあきらかに棄却できる) これでもよい. Pesaran (2007) の表を元に有意水準未満かどうかだけを判断すればよい.
しかし, 全ての問題に対してこのような実装で返された数値が無視できる誤差だとは限らない. IPS 検定の統計量は Dickey Fuller 検定 (Dickey and Fuller 1979) に基づいた方式である. よって検定統計量の正確な計算のためには DF 検定同様にウィーナー過程を含む統計量の計算が必要になるが, これはしょうしょう計算が面倒である. R の他の DF 検定を実装したパッケージも, 実は同様に数表から補間した臨界値とp値を返すものが多い*1.
DF 分布が代数的に計算できず, モンテカルロ法で計算せざるを得ないため, 厳密に一致する値が出せない (乱数による誤差がある) というのが理由の1つで, その点で言えばモンテカルロ法も厳密解でなく近似なので補間であること自体は問題の本質ではなく (DF 分布の数値計算の技術的問題に関しては MacKinnon (2001) などを参考に), 一切注意書きがない plm パッケージのドキュメントに問題があると私は考える. 検定力分析ができないので厳格に計算する必要はないという立場であっても, いずれかが200を超えるケースを全く想定していないのでエラーが発生, は問題だろう*2. 少なくともエラーになる原因が分かりやすいようなエラーメッセージを返すようにしないと不親切だ.
解決法
というわけで, 単純に cipstest 関数の または
が大きすぎるときに発生するインデックスエラーを回避しただけの (正確に計算せずに一番近い値を返す) 修正版をとりあえず用意した. 私のフォークリポジトリからインストールできる.
remotes::install_github("Gedevan-Aleksizde/plm")
または,
git clone git@github.com:Gedevan-Aleksizde/plm.git
install.packages("./", repos = NULL)
この 私家版 plm パッケージをインストールすれば解決するが, 将来本体のほうがバージョンアップしたとか, なんか不安とかいうのであれば R/test_cips.R だけを利用すれば良い. これで従来エラーが出たケースでは警告を発しつつ数表の一番近い値を返すようになる.
ちなみに “time series - How is the augmented Dickey–Fuller test (ADF) table of critical values calculated? - Cross Validatedtime series - How is the augmented Dickey–Fuller test (ADF) table of critical values calculated? - Cross Validated” でも提示されている最もシンプルなタイプの DF 分布は, 以下のように与えられる ( は標準ウィーナー過程). この分布をモンテカルロ法で10回求め, 標準正規分布と比較したものを図 1に掲載する. ここではウィーナー過程を 5,000個の正規乱数で近似した上で 50,000回サンプリングしているがそれでも目に見えて誤差があることがわかる.
シミュレーションとグラフ作成のコードは以下の通り.
require(tidyverse) require(ggthemes) t <- 5000 reps <- 50000 ## it may takes a few minutes.. set.seed(42) d <- list() for(j in 1:10){ d[[j]] <- rep(NA, reps) for (i in 1:reps){ u <- rnorm(t) W <- 1/sqrt(t)*cumsum(u) d[[j]][i] <- (W[t]^2-1)/(2*sqrt(mean(W^2))) } } d <- as_tibble(d, .name_repair = "universal") %>% set_names(1:10) %>% pivot_longer(everything()) %>% mutate(dist = "Dickey Fuller") ggplot(d, aes(x = value, group = name, color = dist, linetype = dist)) + stat_density(size = .5, geom = "line", position = "identity") + stat_function(fun = dnorm, aes(group = "Normal", color = "Normal", linetype = "Normal"), size = 1) + scale_linetype_manual(values = c(3, 1)) + scale_color_pander() + theme_classic(base_size = 14) + theme(legend.title = element_blank(), legend.position = "bottom")
CIPS 検定をモンテカルロで計算するバージョンは気が向いたらやるかもしれないし, やらないかもしれない.
結論 (という名の蛇足とか御意見番とか)
起動時のメッセージにあるように, R は「完全に無保証」なオープンソースであるし, パッケージが CRAN に登録されていることは品質保証に合格したことを表すわけではない*3. そういうのが欲しいなら rOpensci のような品質にコミットしている団体の開発しているものだけを使うか, SAS なり SPSS なり STATA なり昔から使われている有償ソフトウェアを選択すべきだろう (もちろんコミットしているとか昔から使われているという実績とかは絶対の保証ではないが). そういうわけで, 近年 R に限らずデータ分析ソフトウェアを誰でも「自由に」使えるようになっているが, 誰でも「適切に」使えるようになったわけではない. 今回はエラーが出たという分かりやすい理由で不具合を見つけることができたが, プログラムのテストをエラーが出なかったというだけで合格判定するエンジニアはまずいないと思う. よって適切に使用できているかを判断するにはユーザーがパッケージの使い方を知っているだけでなく仕組みを理解することが必要だが, こういう話はともすると「なぜそんなことも知らないのか」「英語のドキュメントくらい読め」「おれはもっと高度な理論を知っている」という本質を外れた知識マウンティングになりがちだ. この文脈での「仕組みを理解する」ということが, 実際にはどこに力点を置くべきなのかはまだまだ議論する余地があるだろう. ということを一昨年くらいから考えている (特に2つ目については knitr 開発者のドキュメントの翻訳といった取り組みを初めている.). 私としても誰でもデータ分析ソフトウェアを本当の意味で自由に使えるようになってほしいとは思っているが, どうやって実現するかは難しい.
参考文献
- Dickey, David A., and Wayne A. Fuller. 1979. “Distribution of the Estimators for Autoregressive Time Series with a Unit Root.” Journal of the American Statistical Association 74 (366a): 427–31. DOI: 10.1080/01621459.1979.10482531.
- MacKinnon, James G. 2001. “Computing Numerical Distribution Functions in Econometrics.” Working Paper 1037. Economics Department, Queen’s University.PDF.
- Pesaran, M. Hashem. 2007. “A Simple Panel Unit Root Test in the Presence of Cross-Section Dependence.” Journal of Applied Econometrics 22 (2): 265–312. doi.org/10.1002/jae.951.
*1:その他のメジャーなパッケージの事情: tseries::adf.test もまた線形補間であるが, ヘルプにその旨を明記している. urca::ur.df も同様. forecast::ndiffs が出力する各種単位根検定は urca の関数を利用して計算しているが, このことは明記されていない. fUnitRoots::adfTest も近似計算であると書いている. R 以外の他のソフトウェア, たとえば stata なんかでも Dickey らの表を補間した値を出力しているとヘルプに書いてあることが多い. また, 近似計算のやりかたにもいろいろあるので, やはりこの仕様は明記すべきだろう.
*2:なお, `plm` には他にもいくつかのパネル単位根検定のバリエーションが用意されているが, 私は使ったことがないためそれらの実装に問題がないかは未確認.
*3:たとえば CRAN 登録の経過を atusy 氏がブログに投稿している: https://blog.atusy.net/2019/06/28/cran-submission/ このように CRAN 側はソフトウェアとしての最低限の動作確認はしているが, 統計手法の実装が適切なのかまでは審査していない
