ARMA項を利用した時の新しい推定手法

短期予測に利用するARMA モデルの推定手法はこれまで条件付き最小二乗法が利用されてきましたが、最新のEViewsには最尤法と一般化最小二乗法が追加され、最尤法がデフォルトになりました。
今回は推定手法の違いがモデル推定に与える影響を考察します。

最新のEViewsでは、モデルにおいてARMA項を利用すると、デフォルトで最尤法による推定を実行するようになりました。
まずはサンプルデータを利用してこれまでのバージョンとの違いを確認してみましょう。サンプルデータは

Program Files/EViews 9/Example Files/EViews9 Manual Data/Chapter22-Time Series Regression

フォルダにあるsowell 22.wf1 を利用します。まだ、最新のEViewsをお持ちでない方は、是非デモ版を入手してください。
このワークファイルにはすでに推定済みの推定式オブジェクトが複数ありますが、それらは無視して新たに推定を行います。

このサンプルデータはSowell (1992a) で利用したデータです。GNP の変化率をARMA(3,2) モデルとして推定します。次の図に示すようにモデルのテキストボックスに入力します。
dlog(gnp) ar(1 to 3) ma(1 to 2) c
もちろん、これは
dlog(gnp) c ar(1) ar(2) ar(3) ma(1) ma(2)
と書いても同じです。定数項のc はdlog(gnp) の右側にあれば、どこにあってもかまいません。

Options タブでMethod 項目がML になっている事を確認します。

OK ボタンをクリックして推定を実行し、推定結果のオブジェクトにeqml という名前を付けます。

Method の項目が ARMA Maximum Likelihood となり、最尤法でパラメータを推定したことが分かります。
収束計算を37回の繰り返し実行したことが分かります。最尤法で推定を行いましたので、残差の分散 σ2 も推定し、推定値の一番下に SIGMASQ として表示しています。
細かい話で恐縮ですが、σ2 は非負ですから、正確に言えば p 値は片側検定になります。つまり、この例で言えば、

show 1-@cnorm(eqml.@tstats(7))

として求めます。
もちろん、この場合、推定画面の表示はt-Statistics になっていますが、正規分布の分布関数を利用します。
実際に計算してみると、ほぼ0であることに変わりはありません。次に同じモデルを GLS(Genelarized Least Squares:一般化最小二乗法) で推定します。
先に示したARMAの項目で GLS を選択します。
さっそく、推定結果を比較してみましょう。オブジェクトの名前は eqgls とします。

こちらには σ2 は存在しません。
推定値を見てみると、定数項とAR項はよく似ているようですが、MA項では有意性や符号が異なっています。
それでは、続けてこれまでデフォルトの手法であった CLS(Condition Least Squares) で推定してみましょう。

推定結果の画面上部を見ると、CLS で推定を行ったが、収束計算に失敗していることが分かります。
ARMAモデルでMA項を利用すると計算が複雑になり、収束計算に失敗してしまうことがあります。
このような結果になってしまったら、残念ですが、この結果をそのまま報告することはできません。
このようなケースでは最適化計算の手法 (Options タブのOptimization の項目) を変更してみましょう。
実際、この場合も Gauss-Newton法を BFGS(Broyden, Fletcher, Goldfarb and Shanno)法に変更すると、収束計算を完了できます。

3つの推定手法

推定手法の概略を説明するために前提となるモデルを決めておきます。
ここでは一般的な ARIMA(p, d, q) モデルとして次のようなモデルを想定します。

1式における L はラグオペレータです。
そしてμt = X′tβ でこれはYt の平均を示します。この時、無条件残差 ut は次のように示すことができます。

そして撹乱項(innovation) は次のような形で示すことにします。

以上の設定で、 ML を選択するとEViewsは無条件残差が正規分布に従うものとして、次に示す正確尤度を最大化します。
詳細は置くとして、ここでは正確尤度の定義だけを示します。

4式の第一項は推定に利用するデータの個数だけで決まります。
第二項のToeplitz 共分散行列Ωの行列式は第三項に含まれるパラメータから求められます。
そこでGLS は次に示す第三項だけを目的関数として、

これを最小化するパラメータを求めます。
従来からEViews がサポートしていた条件付き最小二乗法で推定する場合は次に示す式を対数尤度関数として計算を行うものです。

ここでの計算上のポイントは対数尤度を最大化するために5式の第二項の S (β, ρ, θ) を最小化するところにあります。
何を基準にしてこの3つの推定手法を選択すれば良いでしょうか。
数式を見る限り、パラメータを求める計算に用いる情報量が、 ML, GLS, CLS の順番で徐々に少なくなっている事が分かります。
この事を頭の片隅に置いて、数式を用いた解析的な考察ではなく、乱数発生によるモンテカルロシミュレーションを利用して推定量の分布を考察してみましょう。

数値実験

目的は最尤法(ML)、一般化最小二乗法(GLS)、条件付き最小二乗法(CLS) のうち、どの手法が優れているかを調べることです。
計算速度での評価は行いません。あくまで、設定値に近しい値を推定できるのはどの手法かを見る事で評価します。
最初に変数 Xt とYt の間に次の関係を想定します。

ただし、u0 = e0 = 0 とします。設定値は次のようにします。

以上の設定によって作成したデータから3 つの手法でモデルを推定し、パラメータの分布を考察します。
CLS の項目で解説したように収束計算で失敗することを避けるため、どの推定手法でも最適化にはBFGS を使うことにします。
プログラムの概略は以下の通りです。

プログラム例

tic
wfcreate(wf=unstruct) u 1000
rndseed 12
series xt=nrnd*20
!reps=10000
matrix(!reps,3) arout
matrix(!reps,3) maout
!sigma=1.17
!b0=-1419
!b1=0.85
!rho=0.8
!theta=-0.3
for !i=1 to !reps
smpl @all
series et=nrnd*!sigma
series ut=0
smpl @first @first
ut=et
series yt=!b0+!b1*xt+ut
smpl @first+1 @last

ut=!rho*ut(-1)+et+!theta*et(-1)
yt=!b0+!b1*xt+ut
smpl @all
equation eq01.LS(OPTMETHOD=bfgs,z) yt c xt ar(1) ma(1)
arout(!i,1)=eq01.@coefs(3)
maout(!i,1)=eq01.@coefs(4)
equation eq01.LS(arma=gls,z,OPTMETHOD=bfgs) yt c xt ar(1) ma(1)
arout(!i,2)=eq01.@coefs(3)
maout(!i,2)=eq01.@coefs(4)
equation eq01.LS(arma=cls,z,OPTMETHOD=bfgs) yt c xt ar(1) ma(1)
arout(!i,3)=eq01.@coefs(3)
maout(!i,3)=eq01.@coefs(4)
next
show arout.stats
show maout.stats
toc

プログラム実行結果

結果を見る前にプログラムに関する補足を行います。
この数値実験のサンプルサイズは1000です。
乱数を使って仮想データを作成し、そのデータからパラメータを3つの手法で推定するという処理を1万回行います。
そして1万個の推定値のグラフや記述統計量から分布を考察します。
tic と toc は計算にかかった時間を秒単位でステータスバーに表示するコマンドです。もちろん、プログラムの実行は Quiet モードで行います。
推定コマンドのオプションで z を利用していますが、これは ARMA項のある推定においてバックキャスティングを行わない、というものです。
バックキャスティングについては第三回でご説明します。
この数値実験により得られた1万個の ρ を格納するのが行列オブジェクト arout で、θ は maout に入ります。
最初に ρ の分布を確認します。
ML と GLS はほぼ設定値 0.8 に近くにピークが出現しています。
しかし、CLS は残念ながら平均値は結果として0.71程度ですが、推定量としてみると正規分布とは、かなり異なる分布を示しています。

次にMA項のグラフを示します。この場合も CLS は他の推定手法と比べ、問題のあることが分かります。

最新のEViewsでは推定式においてarma項を利用するとデフォルトで ml による推定が実行されるようになりました。
今回の数値実験を通して、視覚的にそのメリットが体験できたのではないでしょうか。
今回のプログラムを実際に試された方はプログラムを編集し、 β0, β1 の分布についても確認してみましょう。
結果として、 β0, β1 の分布にはほぼ違いはない事が分かります。
次回、第三回では1式における d で整数だけでなく、実数も利用可能な ARFIMAモデルと、今回利用したバックキャスティングという用語の意味とその効果について解説します。

第一回予測平均が一番良い?

EViews 分析機能解説一覧

第三回 ARFIMAモデル