1. 目的

Stataのdfglsコマンドでストアされない臨界値を得る．

2. コード

2.1 データ・検定

ドイツのマクロ経済データより，投資の対数値 ln_inv を検定¹に使用

use http://www.stata-press.com/data/r15/lutkepohl2
dfgls ln_inv

出力：

DF-GLS for ln_inv                                        Number of obs =    80
Maxlag = 11 chosen by Schwert criterion
 
               DF-GLS tau      1% Critical       5% Critical      10% Critical
  [lags]     Test Statistic        Value             Value             Value
------------------------------------------------------------------------------
    11           -2.925           -3.610            -2.763            -2.489
    10           -2.671           -3.610            -2.798            -2.523
    9            -2.766           -3.610            -2.832            -2.555
    8            -3.259           -3.610            -2.865            -2.587
    7            -3.536           -3.610            -2.898            -2.617
    6            -3.115           -3.610            -2.929            -2.646
    5            -3.054           -3.610            -2.958            -2.674
    4            -3.016           -3.610            -2.986            -2.699
    3            -2.071           -3.610            -3.012            -2.723
    2            -1.675           -3.610            -3.035            -2.744
    1            -1.752           -3.610            -3.055            -2.762
 
Opt Lag (Ng-Perron seq t) =  7 with RMSE  .0388771
Min SC   = -6.169137 at lag  4 with RMSE  .0398949
Min MAIC = -6.136371 at lag  1 with RMSE  .0440319

2-2. 臨界値(CV; Critical Value)の取り出し

ラグ次数 = 11の行の[1%CV, 5%CV, 10%CV] = [-3.610, -2.763, -2,489]を取り出すことにする．
1%CVと，5%, 10%CVでdfglsでの計算方法が異なっているため，別々に計算する．

***ここから5%, 10%CVの計算
local N = r(N)
local p = r(maxlag)
tempname rst
mat `rst' = ( -2.550, -20.166, 155.215, 1.133, 9.808, -20.313 \ /*
            */-2.838, -20.328, 124.191, 1.267, 10.530, -24.600 )
local n1 = 1.0/`N'
local pt = `p'/`N'
local cr10 = `rst'[1,1] + `rst'[1,2]*`n1' + `rst'[1,3]*(`n1'^2) + /*
            */ `rst'[1,4]*`pt' + `rst'[1,5]*(`pt'^2) + `rst'[1,6]*(`pt'^3)
local cr05 = `rst'[2,1] + `rst'[2,2]*`n1' + `rst'[2,3]*(`n1'^2) + /*
            */ `rst'[2,4]*`pt' + `rst'[2,5]*(`pt'^2) + `rst'[2,6]*(`pt'^3)
scalar resp10 = `cr10'
scalar resp5 = `cr05'

***ここから1%CVの計算
tempname zt
local N = `N' + `p' + 1
mat `zt' = (-3.77,-3.58,-3.46,-3.48)
if `N' <= 50 {
    local zt1  = `zt'[1,1] 
}
else if `N' <= 100 {
    local zt1  = `zt'[1,1] + (`N'-50)/50 * (`zt'[1,2]-`zt'[1,1])
}
else if `N' <= 200 {
    local zt1  = `zt'[1,2] + (`N'-100)/100 * (`zt'[1,3]-`zt'[1,2])
}
else {
    local zt1  = `zt'[1,4]
}
scalar Zt1  = `zt1'

***ここからCVを格納
mat cv = (Zt1, resp5, resp10)
mat colnames cv = 1%CV 5%CV 10%CV
mat list cv

cv[1,3]
          1%CV        5%CV       10%CV
r1     -3.6104  -2.7633502  -2.4894085

dfgls ln_invの出力結果のラグ次数 = 11の行と一致していることを確認

3. コメント

CV計算詳細はstataマニュアル参照（pdf）．
コードの変数名はdfglsのadoファイルに準ずる．
1%CV計算がCheung & Lai(1995)の表に無かったので少しいい加減になるstata
何度も単位根検定に迫られたときに便利？

dfglsコマンドのデフォルトでは，最大ラグ次数はSchwert (1989)の方法で選択され，トレンド項を含む．いずれもオプションで変更可だが特別変更しない． ↩
観測数，ラグ次数（ここでは11）は，r(N), r(maxlag)にストアされている．Cheung & Lai(1995)の表とは，当該論文のp.413, Table 1(右側)のことで，コードでは行列rstに格納している．当該論文のp.414, (4)式でCVを計算する． ↩
CVはERS(1996)のp.825, Table 1.CのCVを行列ztに格納し，標準的なで補完法で計算される．観測数をN + p + 1としているのは，Cheung & Lai(1995)の方法ではデータセットのサンプルサイズから最大ラグ次数 + 1を差し引いた値が観測数となるが，ERS(1996)の方法ではラグ次数 = 0に基づいているため，データセットのサンプルサイズがそのまま観測数となり，その調整を行っているためである． ↩
実際，dfglsの出力をみると1%CVはラグ次数の値により変化していない． ↩