GARCHモデルのソフトウェア実装とテスト

ホーム > 製品 > nAG数値計算ライブラリ > テクニカルレポート > GARCHモデル

G F LEVY

nAG Ltd, Wilkinson House,
Jordan Hill Road,
Oxford, OX2 8DR, U.K.
(email: george@nag.co.uk)

8 参考文献

概要： 本資料ではFortran 77のGARCH ルーチンのソフトウェア実装について記述しています。ここで考慮されているルーチンは対称GARCHと非対称GARCHの両者を対象とし、ガウス分布と非ガウス分布をもつショックも対象としています。このソフトウェアを使用している様々な例が示されており、モンテカルロシミュレーションによるテスト結果が提示されています。

キーワード：GARCH、ARCH、最尤推定法、ボラティリティ、一般化自己回帰分散不均一、非対称、Fortran 77

1 概要

時間依存性のある分散の数列のモデリングは金融工学の多くの分野で重要です。このテクニカルレポートでは、nAG社で数値計算ライブラリ向けに開発された、一変量の一般化自己回帰条件つき分散不均一（GARCH）ルーチンについて記述しています。ここではFortran 77 ライブラリ [1] の次期バージョン用に開発されたGARCH ソフトウェアに扱いますが、GARCH ソフトウェアのいくつかは既に現バージョンのnAG C ライブラリ [2] に含まれています。ここで記述しているFortran 77 ソフトウェアは任意の値 p と q をもつGARCH(p,q) モデルに使用することができます。GARCH 数列の生成、モデル推定やボラティリティ予測を行うためのルーチンがあります。推定ルーチンはパラメータ推定を返すだけでなく、重要な統計値も返します。たとえば、標準誤差、スコア、計算されたモデルパラメータに対する対数尤度関数のパラメータ値などです。

ソフトウェアの他の機能として以下のものがあります：

回帰－GARCH(p,q) モデル
対称モデル
非対称モデル
ガウス分布と非ガウス分布をもつショック

nAG Fortran 77 と C ライブラリ両者はダイナミックリンクライブラリ（DLL）として実装されているので、Microsoft Windows 環境から簡単に使用することができます。つまり、GARCH ルーチンはExcel、Visual Basic などのMicrosoft ソフトウェアに簡単に組み込むことができます。

2 GARCH モデル

ガウスショックε_t をもつ標準（対称）回帰－GARCH (p,q) モデル [3] [4] [5] は以下の形式で表されます：

このプロセスは、q + 1 個の係数 α_i（i = 0,…,q）、p 個の係数 β_i（i = 1,…,p）、平均 b₀ 、k 個の線形回帰係数b_i（i = 1,…,k）、内生変数 y_t ／外生変数 x_t 、ショックε_t 、条件付き分散 h_t 、そして時間 t までの一連の情報 Ψ_t で記述されます。
p = 0 の場合に GARCH (p,q) モデルはARCH (q) モデルとも呼ばれます。

金融時系列の実証的研究では、金融時系列が負のショック（悪いニュース）の後に条件つき分散 h_t が増加するという特徴があることを示しています。また、ショックの分散はかなりの急尖を示すこともわかりました。標準のガウスGARCH モデルではこれらの効果を得ることができないため、様々な GARCH モデルの拡張版が開発されました [6]。

ここで考慮される対称 GARCH モデルは以下です：

AGARCH (p,q)-type1

AGARCH (p,q)-type2

GJR-GARCH (p,q), または Glosten-Jagannathan-Runkle GARCH [7]

ここでは ε_t > 0 ならば S_t = 1 また ε_t ≥ 0 ならばS_t = 0

EGARCH (p,q), または指数 GARCH

ここでは、と E[| Z_t-i |] は | Z_t-i |の期待値を表します。

これらすべてのモデルでは、ショック ε_t は特定の自由度のガウス分布かスチューデントt-分布のどちらかを持つことができます。

AGARCH-type1 では、追加パラメータ γ により非対称効果がモデル化されます。たとえば、標準 GARCH(1,1) モデルでは、h_t-1が固定の場合、h_t = h (ε_t-1) はε_t-1 = 0 のときに最小となる放物線になります。 ε_t-1 = -γ のときに最小となるよう、追加パラメータ γ を導入すると水平に放物線が移動します。負のショックの後に条件付き分散は γ < 0 を選ぶことにより改善され、ε_t-1 > 0 のときに h (-ε_t-1) > h (ε_t-1) になります。

AGARCH-type2 モデルでは、γ を含めることで、負のショック ε_t-1 の後に h_t を改善することができます。GARCH(1,1) モデルの場合、ε_t-1 > 0 かつ γ < 0 の場合に h (-ε_t-1) > h (ε_t-1) となります。

同様に、GJR-GARCH (1,1) モデルでは、h_t の値はε_t-1 < 0 かつ γ > 0 となるときに増加し、モデルが対称となる場合を上回ります。

EGARCH の場合、項から非対称的な反応が起こります。 EGARCH(1,1) では、 α₁ < 0 の場合、負のショックε_t-1 < 0　は h_t の値を増加させます。これは ln{ h (-Z_t-1) } > ln{h (Z_t-1)}で表わされます。

上記の全ての GARCH プロセスは、パラメータベクトルθで独自に記述されます。この場合、 θ = (b₀,b^T,ω^T), ω^T = (α₀, α₁, α₂,…,α_q,β ₁,β₂,…,β_p,γ) かつ b^T = (b₁,…,b_k)です。

GARCH モデルの実装はここでは全て対数尤度（目的）関数を最大化する θ の値に依存しています。

ここで、T は数列の項数です。

θに対する初期近似から始めて、数値最適化を使用し反復を行って好ましい解を得ることによりこの実装は実現されます。 EGARCH 以外の、全ての GARCH 推定ルーチンでは、目的関数のヤコビアンの解析的表現が最適化の段階で使用されています。

パラメータ推定の標準誤差は既知の結果[8] を使用して計算することができます。その結果とは、θ に対する最尤推定値が、ℑ（フィッシャー情報行列)が以下の式である場合の共分散行列ℑ^-1と平均 θ をもつ漸近的な正規分布となるというものです：

3 GARCH ソフトウェア

ここでは、数列生成、モデルパラメータ推定、回帰GARCH(p,q) 数列の予測に利用できるソフトウェアを記載しています。

3.1 生成

以下のルーチンは、様々な対称及び非対称 GARCH (p,q) 数列から一定数の項を生成します。

仕様

AGARCH-type1: SUBROUTINE G05HKF (DIST, NUM, IP, IQ, THETA, GAMMA, DF, HT, YT, FCALL, RVEC, IFLAG)
AGARCH-type2: SUBROUTINE G05HLF (DIST, NUM, IP, IQ, THETA, GAMMA, DF, HT, YT, FCALL, RVEC, IFLAG)
GJR-GARCH: SUBROUTINE G05HMF (DIST, NUM, IP, IQ, THETA, GAMMA, DF, HT, YT, FCALL, RVEC,IFLAG)
EGARCH: SUBROUTINE G05HNF (DIST, NUM, IP, IQ, THETA, DF, HT, YT, FCALL, RVEC, IFLAG)

ルーチンのパラメータは以下の意味があります：

パラメータ

DIST - CHARACTER*1.
呼び出し時、ε_iに対して使用する分布の種類
DIST = 'N' の場合、正規分布が使用されます。
DIST = 'T' の場合、スチューデントt-分布が使用されます。

NUM - INTEGER
呼び出し時、数列の項数, T

IP - INTEGER.
呼び出し時、移動平均係数の数, p

IQ - INTEGER.
呼び出し時、自己回帰係数の数, q

THETA- DOUBLE PRECISION 配列.
呼び出し時、GARCH モデルのパラメータ

GAMMA - DOUBLE PRECISION.
呼び出し時、GARCH 数列の非対称パラメータ

DF - DOUBLE PRECISION
呼び出し時、スチューデントt-分布の自由度数
DIST = 'N'の場合は、参照されません。

HT - DOUBLE PRECISION 配列.
復帰時、GARCH 数列の条件つき分散, h_i (i = 1,…,T)

YT - DOUBLE PRECISION 配列.
復帰時、GARCH 数列の観測値, ε_i(i= 1,…,T）。

FCALL - LOGICAL.
呼び出し時、FCALL = .TRUE. の場合、新しい数列が生成されます。
そうでない場合は、任意の数列がRVEC の情報を使用して継続されます。

RVEC - DOUBLE PRECISION 配列
呼び出し時、FCALL = .FALSE.の場合、数列を継続するのに必要な情報を配列は含みます。
復帰時、FCALL = .FALSE.の場合、次の呼び出しの中で使用できる情報を配列は含みます。

IFLAG - INTEGER.
エラーインジケータ

3.2 推定

以下のルーチンは様々な対称及び非対称の回帰－GARCH (p,q) 数列のモデルパラメータの推定を行います。

仕様

AGARCH-type1: SUBROUTINE G13FAF (DIST, YT, X, LDX, NUM, IP, IQ, NREG, MN, ISYM, THETA, SE, SC, COVAR, LDC, HP, ETM, HTM, LGF, COPTS, MAXIT, TOL, WORK, LWORK, IFLAG)
AGARCH-type2: SUBROUTINE G13FCF (DIST, YT, X, LDX, NUM, IP, IQ, NREG, MN, THETA, SE, SC, COVAR, LDC, HP, ETM, HTM, LGF, COPTS, MAXIT, TOL, WORK, LWORK, IFLAG)
GJR-GARCH: SUBROUTINE G13FEF (DIST, YT, X, LDX, NUM, IP, IQ, NREG, MN, THETA, SE, SC, COVAR, LDC, HP, ETM, HTM, LGF, COPTS, MAXIT, TOL, WORK, LWORK, IFLAG)
EGARCH: SUBROUTINE G13FGF (DIST, YT, X, LDX, NUM, IP, IQ, NREG, MN, THETA, SE, SC, COVAR, LDC, HP, ETM, HTM, LGF, COPT, MAXIT, TOL, WORK, LWORK, IFLAG)

ルーチンのパラメータは以下の意味があります：

パラメータ

DIST - CHARACTER*1
呼び出し時、ε_iに対して使用する分布の種類。DIST = 'N' の場合、正規分布が使用され、 DIST = 'T' の場合、スチューデントt-分布が使用されます。

YT - DOUBLE PRECISION 配列.
呼び出し時、観測値の数列、ε_i(i, = 1,…,T）。

X - DOUBLE PRECISION 配列.
呼び出し時、X の i 番目の行は時間依存の外生変数ベクトルを含みます。 x_i (i = 1,…,T)　

LDX - INTEGER.
呼び出し時、一次元配列 X

NUM - INTEGER.
呼び出し時、数列の項数, T

IP - INTEGER.
呼び出し時、移動平均係数の数, p

IQ - INTEGER.
呼び出し時、自己回帰係数の数, q

NREG - INTEGER.
呼び出し時、回帰係数の数, k

MN - INTEGER.
呼び出し時、MN = 1 の場合、平均の項b₀ がモデルに含まれます。

ISYM - INTEGER.
呼び出し時、ISYM = 1 の場合、非対称の項γ がモデルに含まれます。（G13FAF にのみ適用）

THETA - DOUBLE PRECISION 配列.
呼び出し時、θ に対する初期パラメータ推定値
復帰時、ベクトル θ に対する推定値 θˆ

SE - DOUBLE PRECISION 配列.
復帰時、θˆ に対する標準誤差

SC - DOUBLE PRECISION 配列.
復帰時、θˆ に対するスコア

COVAR - DOUBLE PRECISION 配列.
復帰時、フィッシャー情報行列の逆の、パラメータ推定θˆ の共分散行列

LDC - INTEGER.
呼び出し時、一次元配列 COVAR

HP - DOUBLE PRECISION
呼び出し時、COPTS(2) = .FALSE. の場合、HP は以前に観測された条件つき分散に使用される値です。

COPTS(1) = .TRUE. の場合、HP は参照されません。
復帰時、COPTS(2) = .TRUE. の場合、HP は以前に観測された条件つき分散の推定値です。

> ET - DOUBLE PRECISION 配列.
復帰時、推定残差, ε_i(i = 1,…,T）

HT - DOUBLE PRECISION 配列.
復帰時、推定条件つき分散, h_i (i = 1,…,T)

LGF - DOUBLE PRECISION
復帰時、 θˆ での尤度関数の値

COPTS - LOGICAL 配列
COPTS(1) = .TRUE. の場合、静的状態になります。そうでない場合は、静的状態にはなりません。
COPTS(2) = .TRUE. の場合、ルーチンは回帰項の初期パラメータ推定を提供します。
そうでない場合は、ユーザが提供します。

MAXIT - INTEGER.
呼び出し時、GARCHパラメータの推定の際に、最適化ルーチンによって使用される最大反復数

TOL - DOUBLE PRECISION
呼び出し時、GARCHパラメータの推定の際に、最適化ルーチンによって使用される公差

WORK - DOUBLE PRECISION 配列, ワークスペース

LWORK - INTEGER.
呼び出し時、作業用配列 WORK のサイズ

IFLAG - INTEGER.
エラーインジケータ

3.3 予測

以下のルーチンは、対称及び非対称 GARCH(p,q) 数列に対するボラティリティの予測を計算します。

仕様

AGARCH-type1: SUBROUTINE G13FBF (NUM, NT, IP, IQ, THETA, GAMMA, CVAR, HT, ET, IFLAG)
AGARCH-type2: SUBROUTINE G13FDF (NUM, NT, IP, IQ, THETA, GAMMA, CVAR, HT, ET, IFLAG)
GJR-GARCH: SUBROUTINE G13FFF (NUM, NT, IP, IQ, THETA, GAMMA, CVAR, HT, ET, IFLAG)
EGARCH: SUBROUTINE G13FHF (NUM, NT, IP, IQ, THETA, CVAR, HT, ET, IFLAG)

ルーチンのパラメータは以下の意味があります：

パラメータ

NUM - INTEGER.
呼び出し時、モデル化された数列からの、配列 HT と ET の項数

NT - INTEGER.
呼び出し時、予測範囲、χ

IP - INTEGER.
呼び出し時、移動平均係数の数, p

IQ - INTEGER.
呼び出し時、自己回帰係数の数, q

THETA - DOUBLE PRECISION 配列.
呼び出し時、GARCH 数列のモデルパラメータ.

GAMMA - DOUBLE PRECISION.
呼び出し時、GARCH 数列の非対称パラメータγ

CVAR - DOUBLE PRECISION 配列.
復帰時,条件つき分散の予測期待値,

HT - DOUBLE PRECISION 配列.
呼び出し時、GARCH (p,q)プロセスの過去の条件つき分散の数列, h_i (i = 1,…,T)

ET - DOUBLE PRECISION 配列.
呼び出し時、GARCH (p,q)プロセスの過去の残差の数列、 ε_i(i = 1,…,T)

IFLAG - INTEGER.
エラーインジケータ

4 使用例

このセクションでは、GARCH ルーチンが実際にどのように使用されるかを示すためにFortran 77 の完全なソースコードを提供しています。それぞれの例は、適切なGARCH推定を用いて任意のGARCH 数列モデルを生成し、ボラティリティ予測を行っています。全ての例で、ガウス分布とスチューデントt-分布のショック両方をもつGARCH数列を考慮しています。

4.1 AGARCH-type1

この例では、以下の二つのモデルが考慮されています：

ガウス分布からのショックと観測値をもつAARCH (3)-type1モデル、y_t = b₀ + ε_t
スチューデントt-分布からのショックと観測値をもつAGARCH (1,2)-type2モデル、

1500個の観測値の数列が、これらのプロセス両方に対して生成されます。そして真値の半分となる初期パラメータ推定を用いてモデル化されます。次に最終的なパラメータ推定が出力され、４段階先のボラティリティ予測が計算されます。

Fortran ソースコード

      
　　　INTEGER NPARMX,NUM
　　　DOUBLE PRECISION ZERO
　　　PARAMETER (NPARMX=10,NUM=1500,ZERO=0.0D0)
　　　INTEGER NUM1,NREGMX,MXNT,NT
　　　PARAMETER (NUM1=3000,NREGMX=10,MXNT=400)
　　　DOUBLE PRECISION FAC1,GAMMA,HP,LGF,MEAN,TOL,XTERM
　　　INTEGER I,IFLAG,IP,IQ,ISYM,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEED
　　　LOGICAL FCALL
　　　DOUBLE PRECISION BX(10),COVAR(NPARMX,NPARMX),ETM(NUM1),
     +                    HT(NUM1+10),HTM(NUM1),PARAM(NPARMX),
     +                    RVEC(40),SC(NPARMX),SE(NPARMX), THETA(NPARMX),
     +                    WK(NUM1*3+NPARMX+NREGMX*NUM1+20*20+1),X(NUM1,10),
     +                    YT(NUM1+10),CVAR(100)
　　　LOGICAL COPTS(2)
　　　CHARACTER*1 DIST
　　　DOUBLE PRECISION DF
　　　EXTERNAL E04UEF,G05HKF,G13FAF,G13FBF,G05CBF
　　　INTRINSIC ABS,DBLE,SIN
　　　WRITE(*,*)'G13FAF Example Program Results'
　　　SEED = 111
　　　NREG = 0
　　　LDX = NUM1
　　　BX(1) = 1.5D0
　　　BX(2) = 2.5D0
　　　BX(3) = 3.0D0
　　　MEAN = 3.0D0
　　　DO 5 I = 1,NUM
　　　　　FAC1 = DBLE(I)*0.01D0
　　　　　X(I,1) = 0.01D0 + 0.7D0*SIN(FAC1)
　　　　　X(I,2) = 0.5D0 + FAC1*0.1D0
　　　　　X(I,3) = 1.0D0
　　5 CONTINUE
　　　ISYM = 1
　　　MN = 1
　　　GAMMA = -0.4D0
　　　IP = 0
　　　IQ = 3
　　　PARAM(1) = 0.8D0
　　　PARAM(2) = 0.6D0
　　　PARAM(3) = 0.2D0
　　　PARAM(4) = 0.1D0
　　　NPAR = 1 + IQ + IP
　　　LWK = NREG*NUM+3*NUM+NPAR+ISYM+MN+NREG+403
　　　FCALL = .TRUE.
　　　IFLAG = 0
　　　DIST = 'N'
　　　CALL G05CBF(SEED)
　　　CALL G05HKF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
     +              FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HKF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
     +              FCALL,RVEC,IFLAG)
　　　IFLAG = -1
　　　DO 10 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 15 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
   15 CONTINUE
         IF (MN.EQ.1) THEN
　　　　　　　YT(I) = MEAN + XTERM + YT(I)
         ELSE
　　　　　　　YT(I) = XTERM + YT(I)
         END IF
   10 CONTINUE
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　COPTS(1) = .TRUE.
　　　COPTS(2) = .TRUE.
　　　MAXIT = 200
　　　TOL = 1.0D-16
　　　DO 12 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
   12 CONTINUE
　　　IF (ISYM.EQ.1) THEN
　　　　　THETA(NPAR+ISYM) = GAMMA*0.5D0
　　　END IF
　　　IFLAG = 0
　　　CALL G13FAF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,ISYM,
     +              THETA,SE,SC,COVAR,
     +              NPARMX,HP,ETM,HTM,LGF,COPTS,MAXIT,TOL,WK,
     +              LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Gaussian distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)
     + 'Parameter estimates Standard errors Correct values'
　　　DO 33 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(I),
     +          SE(I),PARAM(I)
   33 CONTINUE
　　　IF (ISYM.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+1),
     +              SE(NPAR+1),GAMMA
　　　END IF
　　　IF (MN.EQ.1) THEN
           WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM+1),
     +             SE(NPAR+ISYM+1),MEAN
　　　END IF
　　　DO 34 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM+MN+I),
     +      SE(NPAR+ISYM+MN+I),BX(I)
   34 CONTINUE
　　　NT = 4
　　　CALL G13FBF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　WRITE (*,*)
　　　DIST = 'T'
　　　NREG = 2
　　　MN = 1
　　　DF = 4.1D0
　　　IP = 1
　　　IQ = 2
　　　ISYM = 1
　　　GAMMA = -0.2D0
　　　NPAR = IQ + IP + 1
　　　LWK = NREG*NUM+3*NUM+NPAR+ISYM+MN+NREG+404
　　　PARAM(1) = 0.1D0
　　　PARAM(2) = 0.2D0
　　　PARAM(3) = 0.3D0
　　　PARAM(4) = 0.4D0
　　　PARAM(5) = 0.1D0
　　　FCALL = .TRUE.
　　　CALL G05CBF(SEED)
　　　CALL G05HKF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HKF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　CALL G05HKF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　IFLAG = -1
　　　DO 110 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 115 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
  115     CONTINUE
　　　　　IF (MN.EQ.1) THEN
　　　　　　　YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
　　　　　　　YT(I) = XTERM + YT(I)
　　　　　END IF
  110 CONTINUE
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　COPTS(1) = .TRUE.
　　　COPTS(2) = .TRUE.
　　　MAXIT = 200
　　　TOL = 1.0D-16
　　　DO 112 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
  112 CONTINUE
　　　THETA(NPAR+ISYM) = GAMMA*0.5D0
　　　THETA(NPAR+ISYM+1) = DF*0.5D0
　　　CALL G13FAF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,ISYM,
　　 +            THETA,SE,SC,COVAR,NPARMX,HP,ETM,HTM,LGF,
　　 +            COPTS,MAXIT,TOL,WK,LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Student t-distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)
　　 +'Parameter estimates Standard errors Correct values'
　　　DO 133 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(I),
　　 +          SE(I),PARAM(I)
  133 CONTINUE
　　　IF (ISYM.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM),
　　 +          SE(NPAR+ISYM),GAMMA
　　　END IF
　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM+1),
　　 +      SE(NPAR+ISYM+1),DF
　　　IF (MN.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM+1+MN),
　　 +            SE(NPAR+ISYM+1+MN),MEAN
　　　END IF
　　　DO 134 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM+1+MN+I),
　　 +          SE(NPAR+ISYM+1+MN+I),BX(I)
  134 CONTINUE
  199 CONTINUE
　　　NT = 4
　　　CALL G13FBF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　END

出力結果

G13FAF Example Program Results

Gaussian distribution

Parameter estimates  Standard errors  Correct values
　　　　 0.8031            0.0788       0.8000
　　　　 0.6249            0.0570       0.6000
　　　　 0.1803            0.0327       0.2000
　　　　 0.0921            0.0237       0.1000
　　　　-0.5119            0.0682      -0.4000
　　　　 2.9860            0.0324       3.0000

Volatility forecast = 2.8040

Student t-distribution

Parameter estimates  Standard errors  Correct values
         0.0871            0.0230       0.1000
         0.2174            0.0488       0.2000
         0.2736            0.0820       0.3000
         0.3588            0.0788       0.4000
        -0.3240            0.0598      -0.2000
         4.5173            0.5128       4.1000
         3.0182            0.0431       3.0000
         1.4727            0.0265       1.5000
         2.4640            0.0302       2.5000

Volatility forecast = 0.4133

4.2 AGARCH-type2

この例では、以下の二つのモデルが考慮されています：

ガウス分布からのショックと観測値をもつAGARCH (1,1)-type2モデル、
スチューデントt-分布からのショックと観測値をもつAGARCH (1,1)-type2モデル,

1500個の観測値の数列がこれらのプロセス両方に対して生成されます。そして真値の半分となる初期パラメータ推定を用いてモデル化されます。次に最終的なモデルパラメータ推定が出力され、４段階先のボラティリティの予測が計算されます。

Fortran ソースコード

      INTEGER NPARMX,NUM
　　　DOUBLE PRECISION ZERO
　　　PARAMETER (NPARMX=10,NUM=1500,ZERO=0.0D0)
　　　INTEGER NUM1,MXNT,NREGMX,NT
　　　PARAMETER (NUM1=3000,MXNT=400,NREGMX=10)
　　　DOUBLE PRECISION FAC1,GAMMA,HP,LGF,MEAN,TOL,XTERM
　　　INTEGER I,IFLAG,IP,IQ,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEED
　　　LOGICAL FCALL
　　　DOUBLE PRECISION BX(10),COVAR(NPARMX,NPARMX),ETM(NUM1),
　　 +                 HT(NUM1+10),HTM(NUM1),PARAM(NPARMX),
　　 +                 RVEC(40),SC(NPARMX),SE(NPARMX),THETA(NPARMX),
　　 +                 WK(NUM1*3+NPARMX+NREGMX*NUM1+20*20+1),X(NUM1,10),
　　 +                 YT(NUM1+10),CVAR(100)
　　　LOGICAL COPTS(2)
　　　CHARACTER*1 DIST
　　　DOUBLE PRECISION DF
　　　EXTERNAL E04UEF,G05HLF,G13FCF,G13FDF,G05CBF
　　　INTRINSIC ABS,DBLE,SIN
　　　WRITE(*,*)'G13FCF Example Program Results'
　　　SEED = 111
　　　LDX = NUM1
　　　BX(1) = 1.5D0
　　　BX(2) = 2.5D0
　　　BX(3) = 3.0D0
　　　MEAN = 3.0D0
　　　DO 5 I = 1,NUM
　　　　　FAC1 = DBLE(I)*0.01D0
　　　　　X(I,1) = 0.01D0 + 0.7D0*SIN(FAC1)
　　　　　X(I,2) = 0.5D0 + FAC1*0.1D0
　　　　　X(I,3) = 1.0D0
　　5 CONTINUE
　　　MN = 1
　　　NREG = 2
　　　GAMMA = -0.4D0
　　　IP = 1
　　　IQ = 1
　　　NPAR = IQ + IP + 1
　　　LWK = NREG*NUM+3*NUM+NPAR+NREG+MN+404
　　　PARAM(1) = 0.08D0
　　　PARAM(2) = 0.2D0
　　　PARAM(3) = 0.7D0
　　　FCALL = .TRUE.
　　　CALL G05CBF(SEED)
　　　DIST = 'N'
　　　DF = 4.1D0
　　　CALL G05HLF(DIST,300,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HLF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　DO 110 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 120 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
  120     CONTINUE
          IF (MN.EQ.1) THEN
　　　　　　　YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
　　　　　　　YT(I) = XTERM + YT(I)
　　　　　END IF
  110 CONTINUE
　　　IFLAG = -1
　　　DO 130 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
  130 CONTINUE
　　　THETA(NPAR+1) = GAMMA*0.5D0
　　　IF (MN.EQ.1) THEN
           THETA(NPAR+1+MN) = MEAN*0.5D0
　　　END IF
　　　DO 135 I = 1,NREG
　　　　　THETA(NPAR+1+MN+I) = BX(I)*0.5D0
  135 CONTINUE
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　MAXIT = 50
　　　TOL = 1.0D-12
　　　COPTS(1) = .TRUE.
　　　COPTS(2) = .TRUE.
　　　CALL G13FCF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,THETA,
　　 +            SE,SC,COVAR,NPARMX,
　　 +            HP,ETM,HTM,LGF,COPTS,MAXIT,TOL,WK,LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Gaussian distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)'Parameter estimates Standard errors Correct values'
　　　DO 33 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)
   33 CONTINUE
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),GAMMA
　　　IF (MN.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2),
     +          SE(NPAR+2),MEAN
　　　END IF
　　　DO 34 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+MN+1+I),
     +          SE(NPAR+MN+1+I),BX(I)
   34 CONTINUE
　　　NT = 4
　　　CALL G13FDF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　WRITE (*,*)
　　　LWK = NUM1*3 + NPARMX + NREGMX*NUM1 + 1
　　　LDX = NUM1
　　　BX(1) = 1.5D0
　　　BX(2) = 2.5D0
　　　BX(3) = 3.0D0
　　　MEAN = 3.0D0
　　　DO 25 I = 1,NUM
　　　　　FAC1 = DBLE(I)*0.01D0
　　　　　X(I,1) = 0.01D0 + 0.7D0*SIN(FAC1)
　　　　　X(I,2) = 0.5D0 + FAC1*0.1D0
　　　　　X(I,3) = 1.0D0
   25 CONTINUE
　　　MN = 1
　　　NREG = 2
　　　GAMMA = -0.4D0
　　　IP = 1
　　　IQ = 1
　　　NPAR = IQ + IP + 1
　　　LWK = NREG*NUM+3*NUM+NPAR+NREG+MN+405
　　　PARAM(1) = 0.1D0
　　　PARAM(2) = 0.1D0
　　　PARAM(3) = 0.8D0
　　　FCALL = .TRUE.
　　　CALL G05CBF(SEED)
　　　DIST = 'T'
　　　CALL G05HLF(DIST,300,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HLF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
　　 +            FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HLF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,HT,YT,
     +            FCALL,RVEC,IFLAG)
　　　DO 111 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 121 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
  121     CONTINUE
　　　　　IF (MN.EQ.1) THEN
             YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
             YT(I) = XTERM + YT(I)
　　　　　END IF
  111 CONTINUE
　　　IFLAG = -1
　　　DO 131 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
  131 CONTINUE
　　　THETA(NPAR+1) = GAMMA*0.5D0
　　　THETA(NPAR+2) = DF*0.5D0
　　　IF (MN.EQ.1) THEN
　　　　　THETA(NPAR+2+MN) = MEAN*0.5D0
　　　END IF
　　　DO 235 I = 1,NREG
　　　　　THETA(NPAR+MN+2+I) = BX(I)*0.5D0
  235 CONTINUE
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　MAXIT = 100
　　　TOL = 1.0D-12
　　　COPTS(1) = .TRUE.
　　　COPTS(2) = .TRUE.
　　　CALL G13FCF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,THETA,SE,SC,COVAR,
　　 +            NPARMX,HP,ETM,HTM,LGF,COPTS,MAXIT,TOL,WK,LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Student t-distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)'Parameter estimates Standard errors Correct values'
　　　DO 133 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)
  133 CONTINUE
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),GAMMA
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2),SE(NPAR+2),DF
　　　IF (MN.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2+MN),
　　 +          SE(NPAR+2+MN),MEAN
　　　END IF
　　　DO 134 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2+MN+I),
     +          SE(NPAR+2+MN+I),BX(I)
  134 CONTINUE
　　　NT = 4
　　　CALL G13FDF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　END

出力結果

G13FCF Example Program Results

Gaussian distribution

Parameter estimates  Standard errors Correct values
　　　　 0.0835            0.0154      0.0800
　　　　 0.2150            0.0312      0.2000
　　　　 0.6896            0.0324      0.7000
　　　　-0.3757            0.0655     -0.4000
　　　　 3.0453            0.0591      3.0000
　　　　 1.4567            0.0389      1.5000
　　　　 2.4572            0.0445      2.5000

Volatility forecast = 3.0383

Student t-distribution

Parameter estimates  Standard errors  Correct values
　　　　 0.0945            0.0364       0.1000
　　　　 0.0800            0.0264       0.1000
　　　　 0.8197            0.0523       0.8000
　　　　-0.5142            0.1418      -0.4000
　　　　 3.7504            0.3687       4.1000
　　　　 3.0045            0.0631       3.0000
　　　　 1.5321            0.0378       1.5000
　　　　 2.4799            0.0471       2.5000

Volatility forecast = 2.3701

4.3 GJR-GARCH

この例では、以下の二つのモデルが考慮されています：

ガウス分布からのショックと観測値をもつGJR-GARCH (1,1)モデル、
スチューデントt-分布からのショックと観測値をもつGJR-GARCH (1,1)モデル,

2000個の観測値の数列がこれらのプロセス両方に対して生成されます。そして真値の半分となる初期パラメータ推定を用いてモデル化されます。次に最終的なモデルパラメータ推定が出力され、４段階先のボラティリティの予測が計算されます。

Fortran ソースコード

      INTEGER NPARMX,NUM
　　　DOUBLE PRECISION ZERO
　　　PARAMETER (NPARMX=10,NUM=2000,ZERO=0.0D0)
　　　INTEGER NUM1,MXNT,NREGMX,NT
　　　PARAMETER (NUM1=3000,MXNT=400,NREGMX=10)
　　　DOUBLE PRECISION FAC1,GAMMA,HP,LGF,MEAN,TOL,XTERM
　　　INTEGER I,IFLAG,IP,IQ,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEED
　　　LOGICAL FCALL
　　　DOUBLE PRECISION BX(10),COVAR(NPARMX,NPARMX),ETM(NUM1),
　　 +                 HT(NUM1+10),HTM(NUM1),PARAM(NPARMX),
　　 +                 RVEC(40),SC(NPARMX),SE(NPARMX),THETA(NPARMX),
　　 +                 WK(NUM1*3+NPARMX+NREGMX*NUM1+20*20+1),X(NUM1,10),
　　 +                 YT(NUM1+10),CVAR(100)
　　　LOGICAL COPTS(2)
　　　CHARACTER*1 DIST
　　　DOUBLE PRECISION DF
　　　EXTERNAL E04UEF,G05HMF,G13FEF,G13FFF,G05CBF
　　　INTRINSIC ABS,DBLE,SIN
　　　WRITE(*,*)'G13FEF Example Program Results'
　　　SEED = 111
　　　LWK = NUM1*3 + NPARMX + NREGMX*NUM1 + 1
　　　NREG = 0
　　　LDX = NUM1
　　　DF = 5.1D0
　　　GAMMA = 0.1D0
　　　BX(1) = 1.5D0
　　　BX(2) = 2.5D0
　　　BX(3) = 3.0D0
　　　MEAN = 4.0D0
　　　DO 5 I = 1,NUM
　　　　　FAC1 = DBLE(I)*0.01D0
　　　　　X(I,2) = 0.01D0 + 0.7D0*SIN(FAC1)
　　　　　X(I,1) = 0.5D0 + FAC1*0.1D0
　　　　　X(I,3) = 1.0D0
　　5 CONTINUE
　　　MN = 1
　　　NREG = 2
　　　GAMMA = 0.1D0
　　　IP = 1
　　　IQ = 1
　　　NPAR = IQ + IP + 1
　　　PARAM(1) = 0.4D0
　　　PARAM(2) = 0.1D0
　　　PARAM(3) = 0.7D0
　　　FCALL = .TRUE.
　　　DIST = 'N'
　　　CALL G05CBF(SEED)
　　　CALL G05HMF(DIST,200,IP,IQ,PARAM,GAMMA,DF,
　　 +            HT,YT,FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HMF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,
　　 +            HT,YT,FCALL,RVEC,IFLAG)
　　　DO 76 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 77 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
   77     CONTINUE
　　　　　IF (MN.EQ.1) THEN
               YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
　　　　　　　YT(I) = XTERM + YT(I)
　　　　　END IF
   76 CONTINUE
　　　IFLAG = -1
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　COPTS(1) = .TRUE.
　　　COPTS(2) = .TRUE.
　　　MAXIT = 100
　　　TOL = 1.0D-12
　　　DO 81 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
   81 CONTINUE
　　　THETA(NPAR+1) = GAMMA*0.5D0
　　　IF (MN.EQ.1) THEN
           THETA(NPAR+MN+1) = MEAN*0.5D0
　　　END IF
　　　DO 82 I = 1,NREG
　　　　　THETA(NPAR+MN+1+I) = BX(I)*0.5D0
   82 CONTINUE
　　　CALL G13FEF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,THETA,
　　 +            SE,SC,COVAR,NPARMX,
　　 +            HP,ETM,HTM,LGF,COPTS,MAXIT,TOL,WK,LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Gaussian distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)'Parameter estimates Standard errors Correct values'
　　　DO 33 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)
   33 CONTINUE
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),GAMMA
　　　IF (MN.EQ.1) THEN
           WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2),
　　 +           SE(NPAR+2),MEAN
　　　END IF
　　　DO 34 I = 1, NREG
           WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+MN+I+1),
　　 +           SE(NPAR+MN+I+1),BX(I)  
   34 CONTINUE
　　　DIST = 'T'
　　　MEAN = 3.0D0
　　　DO 15 I = 1,NUM
           FAC1 = DBLE(I)*0.01D0
           X(I,2) = 0.01D0 + 0.7D0*SIN(FAC1)
           X(I,1) = 0.5D0 + FAC1*0.1D0
           X(I,3) = 1.0D0
   15 CONTINUE
　　　NT = 4
　　　CALL G13FFF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　WRITE (*,*)
　　　MN = 1
　　　NREG = 2
　　　GAMMA = 0.09D0
　　　IP = 1
　　　IQ = 1
　　　NPAR = IQ + IP + 1
　　　PARAM(1) = 0.05D0
　　　PARAM(2) = 0.1D0
　　　PARAM(3) = 0.8D0
　　　FCALL = .TRUE.
　　　CALL G05CBF(SEED)
　　　CALL G05HMF(DIST,200,IP,IQ,PARAM,GAMMA,DF,
　　 +            HT,YT,FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HMF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,
　　 +            HT,YT,FCALL,RVEC,IFLAG)
　　　CALL G05HMF(DIST,NUM,IP,IQ,PARAM,GAMMA,DF,
　　 +            HT,YT,FCALL,RVEC,IFLAG)
　　　DO 176 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 177 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
  177     CONTINUE
　　　　　IF (MN.EQ.1) THEN
　　　　　　　YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
　　　　　　　YT(I) = XTERM + YT(I)
　　　　　END IF
  176 CONTINUE
　　　IFLAG = -1
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　MAXIT = 100
　　　TOL = 1.0D-14
　　　DO 181 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
  181 CONTINUE
　　　THETA(NPAR+1) = GAMMA*0.5D0
　　　THETA(NPAR+2) = DF*0.5D0
　　　IF (MN.EQ.1) THEN
　　　　　THETA(NPAR+2+MN) = MEAN*0.5D0
　　　END IF
　　　DO 182 I = 1,NREG
　　　　　THETA(NPAR+2+MN+I) = BX(I)*0.5D0
  182 CONTINUE
　　　COPTS(1) = .TRUE.
　　　COPTS(2) = .TRUE.
　　　CALL G13FEF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,THETA,SE,SC,
　　 +            COVAR,NPARMX, HP,ETM,HTM,LGF,COPTS,MAXIT,TOL,WK,LWK,IFLAG)
       WRITE(*,*)
       WRITE(*,*)'Student t-distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)'Parameter estimates Standard errors Correct values'
　　　DO 133 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)
  133 CONTINUE
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),GAMMA
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2),SE(NPAR+2),DF
　　　IF (MN.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2+MN),
　　 +          SE(NPAR+2+MN),MEAN
　　　END IF
　　　DO 134 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2+MN+I),
     +          SE(NPAR+2+MN+I),BX(I) 
  134 CONTINUE
　　　NT = 4
　　　CALL G13FFF(NUM,NT,IP,IQ,THETA,GAMMA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　END

出力結果

G13FEF Example Program Results

Gaussian distribution

Parameter estimates  Standard errors  Correct values
　　　　 0.3706            0.0780       0.4000
　　　　 0.1034            0.0256       0.1000
　　　　 0.7080            0.0413       0.7000
　　　   0.1191            0.0370       0.1000
　　　　 4.0989            0.0950       4.0000
　　　　 1.4255            0.0592       1.5000
　　　　 2.2613            0.0683       2.5000

Volatility forecast = 1.7056

Student t-distribution

Parameter estimates  Standard errors  Correct values
　　　　 0.0377            0.0084       0.0500
　　　　 0.0831            0.0229       0.1000
　　　　 0.8112            0.0260       0.8000
　　　　 0.1161            0.0361       0.0900
　　　　 5.7626            0.6988       5.1000
　　　　 2.9674            0.0363       3.0000
　　　　 1.4891            0.0231       1.5000
　　　　 2.5161            0.0277       2.5000

Volatility forecast = 0.5971

4.4 EGARCH

この例では、以下の二つのモデルが考慮されています。

ガウス分布からのショックと観測値をもつEGARCH (1,1)モデル、
スチューデントt-分布からのショックと観測値をもつEGARCH (1,2)モデル,

Fortran ソースコード

　　　INTEGER NPARMX,NUM
　　　DOUBLE PRECISION ZERO
　　　PARAMETER (NPARMX=10,NUM=1500,ZERO=0.0D0)
　　　INTEGER NUM1,NREGMX,MXNT,NT
　　　PARAMETER (NUM1=3000,NREGMX=10,MXNT=400)
　　　DOUBLE PRECISION FAC1,HP,LGF,MEAN,TOL,XTERM
　　　DOUBLE PRECISION DF
　　　INTEGER I,IFLAG,IP,IQ,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEED
　　　LOGICAL FCALL
　　　CHARACTER*1 DIST
　　　DOUBLE PRECISION BX(10),COVAR(NPARMX,NPARMX),ETM(NUM1),
　　 +                 HT(NUM1+10),HTM(NUM1),PARAM(NPARMX),
　　 +                 RVEC(40),SC(NPARMX),SE(NPARMX),
　　 +                 THETA(NPARMX),
　　 +                 WK(NUM1*3+NPARMX+NREGMX*NUM1+20*20+1),X(NUM1,10),
　　 +                 YT(NUM1+10),CVAR(100)
　　　LOGICAL COPT
　　　EXTERNAL E04UEF,G05HNF,G13FGF,G13FHF,G05CBF
　　　INTRINSIC ABS,DBLE,SIN
　　　WRITE(*,*)'G13FGF Example Program Results'
　　　SEED = 111
　　　LDX = NUM1
　　　BX(1) = 1.5D0
　　　BX(2) = 2.5D0
　　　BX(3) = 3.0D0
　　　MEAN = 3.0D0
　　　DO 5 I = 1,NUM
　　　　　FAC1 = DBLE(I)*0.01D0
　　　　　X(I,1) = 0.01D0 + 0.7D0*SIN(FAC1)
　　　　　X(I,2) = 0.5D0 + FAC1*0.1D0
　　　　　X(I,3) = 1.0D0
　　5 CONTINUE
　　　NREG = 2
　　　MN = 1
　　　IP = 1
　　　IQ = 1
　　　NPAR = IP + 2*IQ + 1
　　　PARAM(1) = 0.1D0
　　　PARAM(2) = -0.3D0
　　　PARAM(3) = 0.1D0
　　　PARAM(4) = 0.9D0
　　　DF = 5.0D0
　　　DIST = 'N'
　　　FCALL = .TRUE.
　　　CALL G05CBF(SEED)
　　　CALL G05HNF(DIST,800,IP,IQ,PARAM,DF,HT,YT,FCALL,
　　 +            RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HNF(DIST,NUM,IP,IQ,PARAM,DF,HT,YT,FCALL,
　　 +            RVEC,IFLAG)
　　　IFLAG = -1
　　　DO 110 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 115 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
  115     CONTINUE
　　　　　IF (MN.EQ.1) THEN
　　　　　　　YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
　　　　　　　YT(I) = XTERM + YT(I)
　　　　　END IF
  110 CONTINUE
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　COPT = .TRUE.
　　　MAXIT = 50
　　　TOL = 1.0D-12
　　　DO 120 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
  120 CONTINUE
　　　IF (MN.EQ.1) THEN
　　　　　THETA(NPAR+MN) = MEAN*0.5D0
　　　END IF
　　　DO 130 I = 1,NREG
　　　　　THETA(NPAR+MN+I) = BX(I)*0.5D0
  130 CONTINUE
　　　LWK = NREG*NUM + 3*NUM + 3
　　　CALL G13FGF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,THETA,
　　 +            SE,SC,COVAR,NPARMX,
　　 +            HP,ETM,HTM,LGF,COPT,MAXIT,TOL,WK,LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Gaussian distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)'Parameter estimates Standard errors Correct values'
　　　WRITE(*,*)
　　　DO 133 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I) 
  133 CONTINUE
　　　IF (MN.EQ.1) THEN
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),MEAN
　　　END IF
　　　DO 134 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+MN+I),SE(NPAR+MN+I),BX(I)
  134 CONTINUE
　　　NT = 4
　　　CALL G13FHF(NUM,NT,IP,IQ,THETA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　WRITE (*,*)
　　　NREG = 2
　　　MN = 1
　　　IP = 1
　　　IQ = 2
　　　NPAR = IP + 2*IQ + 1
　　　PARAM(1) = 0.1D0
　　　PARAM(2) = -0.3D0
　　　PARAM(3) = -0.1D0
　　　PARAM(4) = 0.1D0
　　　PARAM(5) = 0.3D0
　　　PARAM(6) = 0.7D0
　　　DF = 10.0D0
　　　DIST = 'T'
　　　FCALL = .TRUE.
　　　CALL G05CBF(SEED)
　　　CALL G05HNF(DIST,NUM,IP,IQ,PARAM,DF,HT,YT,FCALL,RVEC,IFLAG)
　　　FCALL = .FALSE.
　　　CALL G05HNF(DIST,NUM,IP,IQ,PARAM,DF,HT,YT,FCALL,RVEC,IFLAG)
　　　IFLAG = -1
　　　DO 10 I = 1,NUM
　　　　　XTERM = ZERO
　　　　　DO 15 K = 1,NREG
　　　　　　　XTERM = XTERM + X(I,K)*BX(K)
   15     CONTINUE
　　　　　IF (MN.EQ.1) THEN
　　　　　　　YT(I) = MEAN + XTERM + YT(I)
　　　　　ELSE
　　　　　　　YT(I) = XTERM + YT(I)
　　　　　END IF
   10 CONTINUE
　　　CALL E04UEF('Nolist')
　　　CALL E04UEF('Print Level = 0')
　　　COPT = .TRUE.
　　　MAXIT = 50
　　　TOL = 1.0D-12
　　　DO 20 I = 1,NPAR
　　　　　THETA(I) = PARAM(I)*0.5D0
   20 CONTINUE
　　　THETA(NPAR+1) = DF*0.5D0
　　　IF (MN.EQ.1) THEN
　　　　　THETA(NPAR+1+MN) = MEAN*0.5D0
　　　END IF
　　　DO 30 I = 1,NREG
　　　　　THETA(NPAR+1+MN+I) = BX(I)*0.5D0
   30 CONTINUE
　　　LWK = NREG*NUM + 3*NUM + 3
　　　CALL G13FGF(DIST,YT,X,LDX,NUM,IP,IQ,NREG,MN,THETA,
　　 +            SE,SC,COVAR,NPARMX,
　　 +            HP,ETM,HTM,LGF,COPT,MAXIT,TOL,WK,LWK,IFLAG)
　　　WRITE(*,*)
　　　WRITE(*,*)'Student t-distribution'
　　　WRITE(*,*)
　　　WRITE(*,*)'Parameter estimates Standard errors Correct values'
　　　DO 33 I = 1, NPAR
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)
   33 CONTINUE
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),DF
　　　IF (MN.EQ.1) THEN
　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1+MN),
　　 +      SE(NPAR+1+MN),MEAN
　　　END IF
　　　DO 34 I = 1, NREG
　　　　　WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1+MN+I),
　　 +          SE(NPAR+1+MN+I),BX(I)
   34 CONTINUE
　　　NT = 4
　　　CALL G13FHF(NUM,NT,IP,IQ,THETA,CVAR,HTM,ETM,IFLAG)
　　　WRITE (*,*)
　　　WRITE (*,'(A,F12.4)') 'Volatility forecast = ',CVAR(NT)
　　　END

出力結果

G13FGF Example Program Results

Gaussian distribution

Parameter estimates  Standard errors  Correct values
         0.1153            0.0197       0.1000
        -0.3096            0.0258      -0.3000
         0.1210            0.0382       0.1000
         0.8937            0.0141       0.9000
         2.8624            0.0796       3.0000
         1.4518            0.0445       1.5000
         2.5815            0.0534       2.5000

Volatility forecast = 2.9468

Student t-distribution

Parameter estimates  Standard errors  Correct values
　　　　 0.0992            0.0363       0.1000
　　　  -0.2467            0.0548      -0.3000
　　　　-0.1361            0.0529      -0.1000
　　　　 0.0670            0.0846       0.1000
　　　　 0.3498            0.0680       0.3000
　　　　 0.7298            0.0543       0.7000
　　　　 7.9630            1.4374      10.0000
　　　　 2.9934            0.0734       3.0000
　　　　 1.4935            0.0535       1.5000
　　　　 2.4933            0.0564       2.5000

Volatility forecast = 1.4868

5 ソフトウェアのテスト

GARCHパラメータ推定ソフトウェアが正しく実装されているか確認するためにここで使用されている手法は、既知のGARCH 数列をモデル化するためにGARCHパラメータ推定ソフトウェアを使用するという手法です。それぞれのテストでは、既知の回帰－GARCH (p,q) プロセスは固定のパラメータベクトル θ用に定義されています。モンテカルロシミュレーションは、多くの様々な GARCH 数列を繰り返し生成し、θ推定の平均と分散を計算するために使用されています。それぞれのテストでは、以下が行われています：

初期パラメータ推定値は真値の半分になるようにします
GARCH (1,2) 数列のみモデル化されます

時間依存回帰ベクトル x_t にはすべてのテストにおいて以下の形式をもつ要素があります：

定数要素、 1
線形ランプ要素、
正弦波要素,

GARCH(p,q) 数列のモデリングの難しさは、p とq 両方に依存し、プロセス中にどのくらいボラティリティメモリがあるかにも依存します。β _i パラメータの値が高いと、より多くのボラティリティメモリが生じ、そのため正しくモデル化することが難しくなります。また、モデルパラメータの数が増加すると簡単にモデル化することが難しくなります。なぜなら、数値的最適化を行うための多くの変数があるからです。これは以下の順番で難しいことを示唆しています。ARCH (1), ARCH (2), ARCH (3), GARCH (1,1), GARCH (1,2), GARCH (2,2),…,などです。

GARCH (1,2) モデルはある程度難しいソフトウェアテストであることがわかります。既に述べたとおり、シミュレーションはパラメータ推定 θ の平均と推定標準誤差の平均を計算するのに使用されます。これらの値は、モデルパラメータの真値及び実際の標準誤差と比較されます。

シミュレーションの結果は、セクション6の表1から表20で表されています。「推定値」という列名の最初の列は300 回あるいは200 回のシミュレーションを用いたパラメータ推定の平均を示しています。「推定標準誤差」という列名の２番目の列はGARCH ソフトウェアで計算された標準誤差の平均を示しています。「推定値の標準誤差」という列名の３番目の列はパラメータ推定の実際の標準誤差を示しています。

6 モンテカルロの結果

AGARCH-type1

ここで使用される AGARCH (1,2)-type1 シミュレーションには以下のパラメータがあります：
k = 3
α₀ = 0.2, α₁ = 0.1, α₂ = 0.15, β₁ = 0.7
γ = -0.2
b₀ = 0, b₁ = -1.5, b₂ = 2.5, b₃ = -3.0

系列のショックは、4.10に設定された自由度のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、b_i(i = 1,…,3 )に続いて α₀、α₁,　α₂、β₁、γ　次にdf （もしあれば）という降順に出力されます。

ガウス分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2267	0.1712	0.1147	0.20
0.0873	0.0581	0.0607	0.10
0.1483	0.0903	0.0816	0.15
0.6944	0.1046	0.0793	0.70
-0.2408	0.3024	0.2398	-0.20
-1.4982	0.2586	0.2472	-1.50
2.5562	0.8595	0.7728	2.50
-3.0455	0.6588	0.5912	-3.00

表 1:：400 個の観測値と300 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2022	0.0669	0.0639	0.20
0.0924	0.0411	0.0396	0.10
0.1468	0.0572	0.0526	0.15
0.7050	0.0517	0.0473	0.70
-0.2109	0.1598	0.1380	-0.20
-1.5072	0.0962	0.0971	-1.50
2.5011	0.1619	0.1571	2.50
-3.0001	0.1695	0.1651	-3.00

表 2：1000 個の観測値と 200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2038	0.0467	0.0442	0.20
0.0959	0.0285	0.0283	0.10
0.1526	0.0395	0.0378	0.15
0.6986	0.0357	0.0332	0.70
-0.2076	0.0991	0.0944	-0.20
-1.5004	0.0708	0.0661	-1.50
2.4993	0.0576	0.0559	2.50
-2.9998	0.0952	0.0902	-3.00

表 3：2000 個の観測値と 200 回のシミュレーション

スチューデントt-分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1983	0.0706	0.0676	0.20
0.0937	0.0496	0.0592	0.10
0.1528	0.0684	0.0765	0.15
0.6970	0.0586	0.0582	0.70
-0.2321	0.1954	0.1850	-0.20
4.3246	0.6645	0.6189	4.10
-1.5000	0.0765	0.0812	-1.50
2.4907	0.1270	0.1274	2.00
-2.9950	0.1331	0.1354	-3.00

表 4：1000個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2011	0.0463	0.0456	0.20
0.0975	0.0369	0.0337	0.10
0.1523	0.0485	0.0472	0.15
0.6954	0.0379	0.0394	0.70
-0.2203	0.1338	0.1213	-0.20
4.2278	0.4241	0.4069	4.10
-1.5064	0.0524	0.0517	-1.50
2.4974	0.0451	0.0439	2.00
-2.9987	0.0743	0.0709	-3.00

表 5：2000 個の観測値と 200 回のシミュレーション

AGARCH-type2

ここで使用される AGARCH (1,2)-type2 シミュレーションには以下のパラメータがあります：
k = 3
α₀ = 0.1, α₁ = 0.1, α₂ = 0.15, β₁ = 0.7
γ = -0.1
b₀ = 0, b₁ = -1.5, b₂ = 2.5, b₃ = -3.0

系列のショックは、4.10に設定された自由度のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、b_i(i = 1,…,3 )に続いて α₀、α₁,　α₂、β₁、γ　次にdf （もしあれば）という降順に出力されます。

ガウス分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1238	0.0651	0.0572	0.10
0.1020	0.0636	0.0627	0.10
0.1459	0.0847	0.0856	0.15
0.6778	0.0914	0.0845	0.70
-0.1097	0.1298	0.1116	-0.10
-1.5245	0.1900	0.1763	-1.50
2.4532	0.6022	0.5550	2.50
-2.9564	0.4586	0.4268	-3.00

表 6：400 個の観測値と300 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1081	0.0324	0.0315	0.10
0.1015	0.0436	0.0402	0.10
0.1482	0.0587	0.0538	0.15
0.6918	0.0441	0.0487	0.70
-0.0993	0.0669	0.0613	-0.10
-1.5029	0.0678	0.0675	-1.50
2.4971	0.1083	0.1089	2.50
-2.9935	0.1116	0.1142	-3.00

表 7：1000 個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1048	0.0219	0.0215	0.10
0.1004	0.0310	0.0285	0.10
0.1491	0.0402	0.0379	0.15
0.6960	0.0329	0.0333	0.70
-0.1017	0.0461	0.0425	-0.10
-1.5006	0.0473	0.0459	-1.50
2.4991	0.0356	0.0391	2.50
-2.9996	0.0571	0.0628	-3.00

表 8：2000 個の観測値と200 回のシミュレーション

スチューデントt-分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1103	0.0397	0.0344	0.10
0.0984	0.0548	0.0568	0.10
0.1547	0.0746	0.0752	0.15
0.6856	0.0666	0.0598	0.70
-0.1060	0.0942	0.0848	-0.10
4.2912	0.6639	0.6037	4.10
-1.4998	0.0460	0.0540	-1.50
2.5087	0.0853	0.0882	2.00
-3.0118	0.0860	0.0931	-3.00

表 9：1000 個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1035	0.0236	0.0221	0.10
0.1005	0.0360	0.0347	0.10
0.1528	0.0488	0.0484	0.15
0.6924	0.0417	0.0404	0.70
-0.0998	0.0599	0.0569	-0.10
4.1854	0.4345	0.3988	4.10
-1.4971	0.0324	0.0360	-1.50
2.4986	0.0284	0.0306	2.00
-2.9977	0.0441	0.0494	-3.00

表 10：2000 個の観測値と200 回のシミュレーション

GJR-GARCH

ここで使用される GJR-GARCH (1,2) シミュレーションには以下のパラメータがあります：
k = 3
α₀ = 0.1, α₁ = 0.15, α₂ = 0.2, β₁ = 0.4
γ = 0.1
b₀ = 0, b₁ = -1.5, b₂ = 2.5, b₃ = -3.0

系列のショックは、4.10に設定された自由度数のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、b_i(i = 1,…,3 )に続いて α₀、α₁,　α₂、β₁、γ　次にdf （もしあれば）という降順に出力されます。

ガウス分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1103	0.0410	0.0364	0.10
0.1367	0.0766	0.0774	0.15
0.2003	0.1185	0.1045	0.20
0.3805	0.1435	0.1270	0.40
0.1006	0.0766	0.0759	0.10
-1.4983	0.1010	0.0992	-1.50
2.4982	0.3288	0.3119	2.50
-3.0025	0.2560	0.2397	-3.00

表 11：400 個の観測値と300 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1048	0.0269	0.0221	0.10
0.1451	0.0537	0.0502	0.15
0.1997	0.0752	0.0669	0.20
0.3936	0.0927	0.0784	0.40
0.0958	0.0511	0.0478	0.10
-1.5038	0.0379	0.0387	-1.50
2.4975	0.0645	0.0627	2.50
-2.9970	0.0677	0.0659	-3.00

表 12：1000 個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1018	0.0172	0.0152	0.10
0.1470	0.0387	0.0358	0.15
0.2001	0.0499	0.0476	0.20
0.3972	0.0602	0.0548	0.40
0.1007	0.0385	0.0340	0.10
-1.5016	0.0261	0.0264	-1.50
2.4998	0.0229	0.0224	2.50
-3.0000	0.0363	0.0361	-3.00

表 13：2000 個の観測値と 200 回のシミュレーション

スチューデントt-分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1027	0.0266	0.0285	0.10
0.1392	0.0691	0.0813	0.15
0.2047	0.0932	0.1028	0.20
0.3904	0.0976	0.1131	0.40
0.1158	0.0655	0.1067	0.10
4.2716	0.6603	0.6395	4.10
-1.5013	0.0328	0.0336	-1.50
2.4977	0.0526	0.0564	2.00
-2.9989	0.0549	0.0599	-3.00

表 14：1000 個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.1008	0.0178	0.0170	0.10
0.1469	0.0501	0.0531	0.15
0.2016	0.0716	0.0685	0.20
0.3940	0.0683	0.0663	0.60
0.1067	0.0489	0.0606	0.10
4.2076	0.4360	0.4040	4.10
-1.4982	0.0221	0.0231	-1.50
2.4994	0.0187	0.0188	2.00
-3.0003	0.0299	0.0314	-3.00

表 15：2000 個の観測値と 200 回のシミュレーション

EGARCH

ここで使用される EGARCH (1,2) シミュレーションには以下のパラメータがあります：
k=3, α₀ = 0.3, α₁ = -0.2, α₂ = -0.25,
Φ₁ = 0.1, Φ₂ = 0.15, β₁ = 0.3
b₀ = 0, b₁ = 1.5, b₂ = 2.5, b₃ = 3.0

系列のショックは、4.10に設定された自由度数のガウス分布かスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、b_i(i = 1,…,3 )に続いて α₀、α₁,　α₂、β₁、γ　次にdf （もしあれば）という降順に出力されます。

ガウス分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2786	0.1016	0.0974	0.30
-0.1965	0.0758	0.0738	-0.20
-0.2481	0.0958	0.0909	-0.25
0.0699	0.1186	0.1228	0.10
0.1346	0.1445	0.1285	0.15
0.3045	0.2213	0.1841	0.30
1.5181	0.1855	0.1928	1.50
2.5622	0.6413	0.6017	2.50
2.9430	0.4929	0.4626	3.00

表 16：400 個の観測値と300 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2958	0.0709	0.0665	0.30
-0.1954	0.0471	0.0462	-0.20
-0.2533	0.0595	0.0571	-0.25
0.0848	0.0768	0.0789	0.15
0.1432	0.0850	0.0812	0.10
0.2914	0.1437	0.1247	0.30
1.4970	0.0747	0.0756	1.50
2.4969	0.1245	0.1231	2.50
3.0004	0.1301	0.1309	3.00

表 17：1000 個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2982	0.0452	0.0471	0.30
-0.1980	0.0317	0.0326	-0.20
-0.2540	0.0381	0.0401	-0.25
0.0932	0.0508	0.0533	0.15
0.1491	0.0591	0.0545	0.10
0.2967	0.0895	0.0870	0.30
1.5017	0.0528	0.0515	1.50
2.4999	0.0442	0.0431	2.50
2.9984	0.0748	0.0702	3.00

表 18：2000 個の観測値と200 回のシミュレーション

スチューデントt-分布

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2811	0.1022	0.0924	0.30
-0.2071	0.0628	0.0639	-0.20
-0.2452	0.0786	0.0794	-0.25
0.0865	0.0929	0.0903	0.10
0.1464	0.0932	0.0924	0.15
0.3223	0.1738	0.1526	0.30
4.3119	0.6819	0.4159	4.10
1.4983	0.0668	0.0738	1.50
2.4871	0.1142	0.1210	2.50
3.0094	0.1195	0.1267	3.00

表 19：1000 個の観測値と200 回のシミュレーション

推定値	推定標準誤差	推定値の標準誤差	正しい値
0.2859	0.0702	0.0811	0.30
-0.2016	0.0426	0.0493	-0.20
-0.2483	0.0594	0.0706	-0.25
0.0926	0.0652	0.0657	0.10
0.1512	0.0633	0.0801	0.15
0.3154	0.1174	0.1283	0.30
4.2281	0.4348	0.3854	4.10
1.4942	0.0457	0.0539	1.50
2.4976	0.0413	0.0451	2.50
3.0001	0.0681	0.0734	3.00

表 20：2000 個の観測値と200 回のシミュレーション

結論

表1 から表20 で表されるシミュレーションの結果は、GARCH 生成及び推定 Fortran 77 ソフトウェアが予想通り動作していることを示しています。（予測ソフトウェアは別の場所でテストされています。[9] ）信頼できるパラメータ推定と関連する標準誤差は少なくとも300の観測値がGARCHモデルに含まれるときにのみ得られることがわかりました。より難しいGARCHモデル（より多くの推定するパラメータをもつモデルや高いボラティリティメモリをもつモデル）が一貫性のあるパラメータ推定を得るにはさらに多くの観測値が必要であることがわかりました。

このレポートではMark 20 の nAG Fortran 77 ライブラリ用に開発されたGARCH ソフトウェアの現在の状況を述べています。今後変更／改善が行われるため、nAG数値計算ライブラリにいずれは含まれるGARCHソフトウェアへのあくまでも手引きとしてとらえることができます。今後行われる可能性のあるソフトウェアへの改善には以下が含まれます：

モデルパラメータのいくつかを修正できるように（つまり、最大尤度では推定を行わない）GARCH 推定ルーチンの変更
一般化誤差分布などの他の非ガウスショック
GARCH-M などの他の一変量GARCH モデル
多変量GARCH モデル
回帰-GARCH (p,q) からARMA (p1,q1)-GARCH (p2,q2)へのGARCHモデルの一般化

言語／ユーザインターフェースに関して行われる可能性のある変更には以下があります：

ユーザインターフェースを改善するためのメモリ割り当てなどの近代的な Fortran の特徴を含めること。Fortran 77 はメモリ割り当てをサポートしていないため、現行のGARCH ソフトウェアではどのくらいワークスペースが必要かをユーザが計算する必要があります。またルーチンは20個の推定パラメータしか許されていません；これでほとんどの要件を満たさなければなりません。内部のメモリ割り当てを使用しているnAG C GARCH ソフトウェアはこれらの制限がありませんので、簡単に使用できます。
メモリ割り当てを実行し任意パラメータなどをもつVisual Basic ラッパーの作成。例えば、これらは対話形式での利用という利点があり、ユーザを生のFortran 77コードから保護するGARCH Microsoft Excelアドインの形式をとっています。
全てのGARCH 関数を含む C++ クラスの作成。これは任意／デフォルトの関数パラメータをもち、ワークスペースを隠蔽します。

結論として、nAG 数値計算ライブラリのPC やワークステーションの実装により、時系列パッケージからの該当する関数の呼び出しよりもGARCH 関数の呼び出しのほうがソフトウェア開発者にとって簡単に行うことができると言えます。時系列パッケージは対話型インターフェースの利点がありますが、大規模システムの構成要素として時系列パッケージをプログラムで呼び出すのは容易ではありません（あるいは計算効率が良くはありません）。現在のnAGの提供形態は新規のアプリケーションに個別のGARCHルーチンを組み込みたいと考えるソフトウェア開発者によく適しています。

7 謝辞

ご支援、ご助言をいただきました Geoff Morgan、そして原稿を校正いただきましたAnneTrefethen と Neil Swindells に心より感謝申し上げます。

全ての商標は商標登録されています。

8 参考文献

[1] nAG Ltd, The Fortran 77 Library Mark 20, nAG Ltd, Oxford, in preparation
[2] nAG Ltd, The nAG C Library Mark 6, nAG Ltd, Oxford, 2000.
[3] Engle R, Autoregressive Conditional Heteroskedasticity with Estimates of the Variance of United Kingdom Inflation, Econometrica, 50, 987-1008, 1982
[4] Bollerslev T, Generalised Autoregressive Conditional Heteroskedasticity, Journal of Econometrics, 31, 307-27, 1986
[5] Hamilton J, Time Series Analysis, Princeton University Press, 1994.
[6] Engle R, and Ng V, Measuring and Testing the Impact of News on Volatility, Journal of Finance, 48, 1749-1777, 1993
[7] Glosten L, Jagannathan R, and Runkle D, Relationship between the Expected Value and the Volatility of the Nominal Excess Return on Stocks, Journal of Finance, 48, 1779-1801, 1993.
[8] Silvey S D, Statistical Inference, Monographs on Applied Probability and Statistics, Chapman and Hall, 1975
[9] nAG Ltd, The Fortran 77 Library Mark 20 GARCH stringent test programs, nAG Ltd, Oxford, 2000

GARCHモデルの : ソフトウェア実装とテスト

テクニカルレポート

G F LEVY

nAG Ltd, Wilkinson House,
Jordan Hill Road,
Oxford, OX2 8DR, U.K.
(email: george@nag.co.uk)

目次

1 概要

2 GARCHモデル

3 GARCHソフトウェア

3.1 生成

3.2 推定

3.3 予測

4 使用例

4.1 AGARCH-type1

4.2 AGARCH-type2

4.3 GJR-GARCH

4.4 EGARCH

5 ソフトウェアのテスト

6 モンテカルロの結果

7 謝辞

8 参考文献

1 概要

2 GARCH モデル

3 GARCH ソフトウェア

3.1 生成

3.2 推定

3.3 予測

4 使用例

4.1 AGARCH-type1

4.2 AGARCH-type2

4.3 GJR-GARCH

4.4 EGARCH

5 ソフトウェアのテスト

6 モンテカルロの結果

結論

7 謝辞

8 参考文献

GARCHモデルの : ソフトウェア実装とテスト

テクニカルレポート

G F LEVY nAG Ltd, Wilkinson House, Jordan Hill Road, Oxford, OX2 8DR, U.K. (email: george@nag.co.uk)

目次

1 概要

2 GARCH モデル

3 GARCH ソフトウェア

3.1 生成

3.2 推定

3.3 予測

4 使用例

4.1 AGARCH-type1

4.2 AGARCH-type2

4.3 GJR-GARCH

4.4 EGARCH

5 ソフトウェアのテスト

6 モンテカルロの結果

結論

7 謝辞

8 参考文献

G F LEVY

nAG Ltd, Wilkinson House,
Jordan Hill Road,
Oxford, OX2 8DR, U.K.
(email: george@nag.co.uk)