Stataで単位根検定(DF-GLS)の臨界値を得る

#1. 目的
 - Stataのdfglsコマンドでストアされない臨界値を得る．

#2. コード
##2.1 データ・検定
 - ドイツのマクロ経済データより，投資の対数値 ln_inv を検定[^1]に使用

```
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
```
 - 5%, 10%CV計算は，[Cheung & Lai(1995)][1]による．
  - 観測数，ラグ次数，[Cheung & Lai(1995)][1]の表を用いる[^2]．
 - 1%CV計算[^3]は，[ERS(1996)][2]による．
  - なぜなら，**[Cheung & Lai(1995)][1]の表に計算に必要な情報が無いため**
  - [ERS(1996)][2]の方法はラグ次数 = 0を前提にCV計算されており，[Cheung & Lai(1995)][1]はラグ次数によりCVが影響をうけると主張している．したがって，dfglsが出力する1%CVは[ERS(1996)][2]による[Cheung & Lai(1995)][1]の値の代理である[^4]．
 - 出力：

```
cv[1,3]
          1%CV        5%CV       10%CV
r1     -3.6104  -2.7633502  -2.4894085
```
  - dfgls ln_invの出力結果のラグ次数 = 11の行と一致していることを確認

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

[1]:https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1468-0084.1995.mp57003008.x
[2]:https://www.jstor.org/stable/2171846?seq=1#page_scan_tab_contents
[3]:https://www.stata.com/manuals/tsdfgls.pdf
[4]:https://amstat.tandfonline.com/doi/abs/10.1080/07350015.1989.10509723#.XGV-6Fz7TDc

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