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

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

1
TR4/2000

NP3533

概要： 本資料では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 モデル

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

されます：
,

∑

|

(0, )

∑

2
TR4/2000

NP3533

このプロセスは、q + 1 個の係数

（i =0,…,q）、p 個の係数

（i 1,…,k）、内生変数の線形回帰係数そして時間 t までの一連の情報／外生変数（i 1,…, p）、平均

、ショック

、k 個

、条件付き分散、で記述されます。

p = 0 の場合に GARCH(p,q) モデルはARCH(q) モデルとも呼ばれます。

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

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

AGARCH(p,q)-type1∑∑

AGARCH(p,q)-type2

∑

|

|

∑

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

∑
ここでは

∑

0 ならば

1 また0 ならば

0

EGARCH(p,q), または指数 GARCH
∑
ここでは、

∑

√

と E|||| は |

これらすべてのモデルでは、ショック

[|

|]

∑

(

)

|の期待値を表します。
は特定の自由度のガウス分布かスチューデントt-分布のどちらかを持つことができます。

3
TR4/2000

NP3533

AGARCH-type1 では、追加パラメータ γ により非対称効果がモデル化されます。たとえば、
標準 GARCH(1,1) モデルでは、が固定の場合、はなる放物線になります。

のときに最小となるよう、追加パラメータ γ を導入すると水平に放物線が移動します。負のショックの後に条件付き分散は

され、

のときに

= 0のときに最小とを選ぶことにより改善になります。

AGARCH-type2 モデルでは、γ を含めることで、負のショック

かつとができます。GARCH(1,1) モデルの場合、の後にを改善するこの場合にとなります。

同様に、GJR-GARCH(1,1) モデルでは、

の値は

かつ

となるときに増加し、モデルが対称となる場合を上回ります。
EGARCH の場合、項 ∑から非対称的な反応が起こります。 EGARCH(1,1) では、の場合、負のショック ln

はの値を増加させます。これは lnで表わされます。

上記の全ての GARCH プロセスは、パラメータベクトル で独自に記述されます。この場合、

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

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

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

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

[∑

∂
∂ ∂

]

4
TR4/2000

NP3533

3GARCH ソフトウェア

ここでは、数列生成、モデルパラメータ推定、回帰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)EGARCHSUBROUTINE G05HNF (DIST, NUM, IP, IQ, THETA, DF, HT, YT, FCALL, RVEC, IFLAG)

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

パラメータ
DIST - CHARACTER*1.
呼び出し時、に対して使用する分布の種類

DIST = 'N' の場合、正規分布が使用されます。
DIST = 'T' の場合、スチューデントt-分布が使用されます。
NUM - INTEGER
呼び出し時、数列の項数, T
IP - INTEGER.
呼び出し時、移動平均係数の数, p
IQ - INTEGER.
呼び出し時、自己回帰係数の数, q
THETA- DOUBLE PRECISION 配列.
呼び出し時、GARCH モデルのパラメータ
GAMMA - DOUBLE PRECISION.
呼び出し時、GARCH 数列の非対称パラメータ5
TR4/2000

NP3533

DF - DOUBLE PRECISION
呼び出し時、スチューデントt-分布の自由度数
DIST = 'N'の場合は、参照されません。
HT - DOUBLE PRECISION 配列.
復帰時、GARCH 数列の条件つき分散,(i 1,…,T)

YT - DOUBLE PRECISION 配列.
復帰時、GARCH 数列の観測値,(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)EGARCHSUBROUTINE G13FGF (DIST, YT, X, LDX, NUM, IP, IQ, NREG, MN, THETA, SE, SC, COVAR,
LDC, HP, ETM, HTM,

6
TR4/2000

NP3533

LGF, COPT, MAXIT, TOL, WORK, LWORK, IFLAG)
ルーチンのパラメータは以下の意味があります：
パラメータ
DIST - CHARACTER*1
呼び出し時、に対して使用する分布の種類。DIST = 'N' の場合、正規分布が使用され、DIST = 'T' の場合、スチューデントt-分布が使用されます。
YT - DOUBLE PRECISION 配列.
呼び出し時、観測値の数列、(i 1,…, T）。

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

(i 1,…,T)

LDX - INTEGER.
呼び出し時、一次元配列 X
NUM - INTEGER.
呼び出し時、数列の項数, T
IP - INTEGER.
呼び出し時、移動平均係数の数, p
IQ - INTEGER.
呼び出し時、自己回帰係数の数, q
NREG - INTEGER.
呼び出し時、回帰係数の数, k
MN - INTEGER.
呼び出し時、MN = 1 の場合、平均の項がモデルに含まれます。

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

THETA - DOUBLE PRECISION 配列.
呼び出し時、に対する初期パラメータ推定値

復帰時、ベクトルに対する推定値 ̂

SE - DOUBLE PRECISION 配列.
復帰時、 ̂ に対する標準誤差
SC - DOUBLE PRECISION 配列.
復帰時、 ̂ に対するスコア
COVAR - DOUBLE PRECISION 配列.
復帰時、フィッシャー情報行列の逆の、パラメータ推定 ̂ の共分散行列
LDC - INTEGER.
呼び出し時、一次元配列 COVAR
HP - DOUBLE PRECISION

7
TR4/2000

NP3533

呼び出し時、COPTS(2) = .FALSE. の場合、HP は以前に観測された条件つき分散に使用される値です。
COPTS(1) = .TRUE. の場合、HP は参照されません。
復帰時、COPTS(2) = .TRUE. の場合、HP は以前に観測された条件つき分散の推定値です。
ET - DOUBLE PRECISION 配列.
復帰時、推定残差,(i 1,…, T）

HT - DOUBLE PRECISION 配列.
復帰時、推定条件つき分散,(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)EGARCHSUBROUTINE G13FHF (NUM, NT, IP, IQ, THETA, CVAR, HT, ET, IFLAG)

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

8
TR4/2000

NP3533

パラメータ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)プロセスの過去の条件つき分散の数列,(i 1,…,T)

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

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

9
TR4/2000

NP3533

使用例

4

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

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

ガウス分布からのショックと観測値をもつAARCH(3)-type1モデル、



スチューデントt-分布からのショックと観測値をもつAGARCH(1,2)-type2モデル、1500個の観測値の数列が、これらのプロセス両方に対して生成されます。そして真値の半分となる初期パラメータ推定を用いてモデル化されます。 次に最終的なパラメータ推定が出力され、４段階先のボラティリティ予測が計算されます。

Fortran ソースコード
INTEGER NPARMX,NUM
DOUBLE PRECISION ZEROPARAMETER (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,XTERMINTEGER I,IFLAG,IP,IQ,ISYM,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEEDLOGICAL FCALLDOUBLE 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 DISTDOUBLE 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

10
TR4/2000

NP3533

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.1D0X(I,3) = 1.0D0
5 CONTINUEISYM = 1
MN = 1GAMMA = -0.4D0IP = 0
IQ = 3
PARAM(1) = 0.8D0PARAM(2) = 0.6D0PARAM(3) = 0.2D0PARAM(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)ELSEYT(I) = XTERM + YT(I)END IF
10 CONTINUECALL E04UEF('Nolist')
CALL E04UEF('Print Level = 0')
COPTS(1) = .TRUE.
COPTS(2) = .TRUE.
MAXIT = 200

11
TR4/2000

NP3533

TOL = 1.0D-16
DO 12 I = 1,NPAR
THETA(I) = PARAM(I)*0.5D0
12 CONTINUEIF (ISYM.EQ.1) THEN
THETA(NPAR+ISYM) = GAMMA*0.5D0
END IFIFLAG = 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, NPARWRITE(*,'(F16.4,F18.4,F13.4)') THETA(I),+SE(I),PARAM(I)33 CONTINUEIF (ISYM.EQ.1) THEN
WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+1),+

SE(NPAR+1),GAMMAEND IF
IF (MN.EQ.1) THEN
WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM+1),+

SE(NPAR+ISYM+1),MEAN

END IFDO 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 = 2ISYM = 1
GAMMA = -0.2D012
TR4/2000

NP3533

NPAR = IQ + IP + 1
LWK = NREG*NUM+3*NUM+NPAR+ISYM+MN+NREG+404
PARAM(1) = 0.1D0
PARAM(2) = 0.2D0PARAM(3) = 0.3D0PARAM(4) = 0.4D0PARAM(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)ELSEYT(I) = XTERM + YT(I)END IF

110 CONTINUECALL 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.5D0112 CONTINUE
THETA(NPAR+ISYM) = GAMMA*0.5D0
THETA(NPAR+ISYM+1) = DF*0.5D0CALL 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(*,*)

13
TR4/2000

NP3533

+ 'Parameter estimates Standard errors Correct values'
DO 133 I = 1, NPARWRITE(*,'(F16.4,F18.4,F13.4)') THETA(I),+SE(I),PARAM(I)133 CONTINUEIF (ISYM.EQ.1) THEN
WRITE(*,'(F16.4,F18.4,F13.4)') THETA(NPAR+ISYM),+

SE(NPAR+ISYM),GAMMAEND 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 IFDO 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.80310.62490.18030.0921-0.51192.98600.07880.05700.03270.02370.06820.03240.80000.60000.20000.1000-0.4000
3.0000Volatility forecast = 2.8040
Student t-distribution
Parameter estimates Standard errors Correct values
0.08710.21740.2736
0.3588-0.3240
4.51733.01820.02300.04880.08200.07880.05980.51280.04310.10000.20000.30000.4000-0.2000
4.10003.000014
TR4/2000

NP3533

1.47272.46400.02650.03021.50002.5000Volatility forecast = 0.4133

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



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



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

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

Fortran ソースコード
INTEGER NPARMX,NUM
DOUBLE PRECISION ZEROPARAMETER (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,XTERMINTEGER I,IFLAG,IP,IQ,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEEDLOGICAL FCALLDOUBLE 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 DISTDOUBLE 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)

15
TR4/2000

NP3533

X(I,2) = 0.5D0 + FAC1*0.1D0X(I,3) = 1.0D0
5 CONTINUE
MN = 1
NREG = 2
GAMMA = -0.4D0IP = 1
IQ = 1
NPAR = IQ + IP + 1
LWK = NREG*NUM+3*NUM+NPAR+NREG+MN+404
PARAM(1) = 0.08D0
PARAM(2) = 0.2D0PARAM(3) = 0.7D0FCALL = .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)ELSEYT(I) = XTERM + YT(I)END IF

110 CONTINUEIFLAG = -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 IFDO 135 I = 1,NREG
THETA(NPAR+1+MN+I) = BX(I)*0.5D0
135 CONTINUECALL E04UEF('Nolist')
CALL E04UEF('Print Level = 0')

16
TR4/2000

NP3533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, NPARWRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)33 CONTINUEWRITE(*,'(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 IFDO 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.1D0X(I,3) = 1.0D0
25 CONTINUE
MN = 1
NREG = 2
GAMMA = -0.4D0IP = 1
IQ = 1
NPAR = IQ + IP + 1

17
TR4/2000

NP3533

LWK = NREG*NUM+3*NUM+NPAR+NREG+MN+405
PARAM(1) = 0.1D0
PARAM(2) = 0.1D0
PARAM(3) = 0.8D0FCALL = .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)ELSEYT(I) = XTERM + YT(I)END IF

111 CONTINUEIFLAG = -1
DO 131 I = 1,NPAR
THETA(I) = PARAM(I)*0.5D0131 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 IFDO 235 I = 1,NREG
THETA(NPAR+MN+2+I) = BX(I)*0.5D0235 CONTINUECALL 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)18
TR4/2000

NP3533

WRITE(*,*)
WRITE(*,*)'Student t-distribution'
WRITE(*,*)
WRITE(*,*)'Parameter estimates Standard errors Correct values'
DO 133 I = 1, NPARWRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)133 CONTINUEWRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),GAMMAWRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2),SE(NPAR+2),DFIF (MN.EQ.1) THEN
WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2+MN),+

SE(NPAR+2+MN),MEAN

END IFDO 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.08350.21500.6896-0.3757
3.04531.45672.45720.01540.03120.03240.06550.05910.03890.04450.08000.20000.7000-0.4000
3.00001.50002.5000Volatility forecast = 3.0383
Student t-distribution
Parameter estimates Standard errors Correct values
0.09450.08000.8197-0.51423.75043.00451.53212.47990.03640.02640.05230.14180.36870.06310.03780.04710.10000.10000.8000-0.4000
4.10003.00001.50002.5000Volatility forecast = 2.3701

19
TR4/2000

NP3533

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



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



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

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

Fortran ソースコード
INTEGER NPARMX,NUM
DOUBLE PRECISION ZEROPARAMETER (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,XTERMINTEGER I,IFLAG,IP,IQ,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEEDLOGICAL FCALLDOUBLE 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 DISTDOUBLE 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

20
TR4/2000

NP3533

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.1D0X(I,3) = 1.0D0
5 CONTINUE
MN = 1
NREG = 2
GAMMA = 0.1D0
IP = 1
IQ = 1
NPAR = IQ + IP + 1
PARAM(1) = 0.4D0PARAM(2) = 0.1D0
PARAM(3) = 0.7D0FCALL = .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)ELSEYT(I) = XTERM + YT(I)END IF76 CONTINUEIFLAG = -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.5D081 CONTINUE
THETA(NPAR+1) = GAMMA*0.5D021
TR4/2000

NP3533

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.5D082 CONTINUECALL 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, NPARWRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)33 CONTINUEWRITE(*,'(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 IFDO 34 I = 1, NREG
WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+MN+I+1),+

SE(NPAR+MN+I+1),BX(I)34 CONTINUEDIST = '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.1D0X(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.1D022
TR4/2000

NP3533

PARAM(3) = 0.8D0FCALL = .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)ELSEYT(I) = XTERM + YT(I)END IF

176 CONTINUEIFLAG = -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 IFDO 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, NPAR23
TR4/2000

NP3533

WRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)133 CONTINUEWRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+1),SE(NPAR+1),GAMMAWRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2),SE(NPAR+2),DFIF (MN.EQ.1) THEN
WRITE(*,'(F16.4,F18.4,F16.4)') THETA(NPAR+2+MN),+

SE(NPAR+2+MN),MEAN
END IFDO 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.37060.10340.70800.1191
4.09891.42552.26130.07800.0256
0.04130.03700.09500.05920.06830.40000.10000.70000.10004.00001.50002.5000Volatility forecast = 1.7056
Student t-distribution
Parameter estimates Standard errors Correct values
0.03770.08310.81120.1161
5.76262.96741.48912.51610.00840.02290.02600.0361
0.69880.03630.02310.02770.05000.10000.80000.09005.10003.00001.50002.5000Volatility forecast = 0.5971

24
TR4/2000

NP3533

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



ガウス分布からのショックと観測値をもつEGARCH(1,1)モデル、



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

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

Fortran ソースコード
INTEGER NPARMX,NUM
DOUBLE PRECISION ZEROPARAMETER (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,XTERMDOUBLE PRECISION DF
INTEGER I,IFLAG,IP,IQ,K,LDX,LWK,MAXIT,MN,NPAR,NREG,SEEDLOGICAL FCALL
CHARACTER*1 DISTDOUBLE 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.1D025
TR4/2000

NP3533

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.3D0PARAM(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)ELSEYT(I) = XTERM + YT(I)END IF

110 CONTINUECALL 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.5D0120 CONTINUE
IF (MN.EQ.1) THEN
THETA(NPAR+MN) = MEAN*0.5D0
END IFDO 130 I = 1,NREG
THETA(NPAR+MN+I) = BX(I)*0.5D0
130 CONTINUE

26
TR4/2000

NP3533

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, NPARWRITE(*,'(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 IFDO 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.3D0PARAM(3) = -0.1D0
PARAM(4) = 0.1D0
PARAM(5) = 0.3D0PARAM(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)27
TR4/2000

NP3533

15

CONTINUE
IF (MN.EQ.1) THEN
YT(I) = MEAN + XTERM + YT(I)ELSEYT(I) = XTERM + YT(I)END IF

10 CONTINUECALL 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.5D030 CONTINUELWK = 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, NPARWRITE(*,'(F16.4,F18.4,F16.4)') THETA(I),SE(I),PARAM(I)33 CONTINUEWRITE(*,'(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 IFDO 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)28
TR4/2000

NP3533

END

出力結果
G13FGF Example Program Results
Gaussian distribution
Parameter estimates Standard errors Correct values
0.11530.01970.1000-0.30960.0258-0.3000
0.12100.03820.10000.89370.01410.90002.86240.07963.00001.45180.04451.50002.58150.05342.5000Volatility forecast = 2.9468
Student t-distribution
Parameter estimates Standard errors Correct values
0.0992-0.2467-0.1361
0.06700.34980.72987.96302.9934
1.49352.49330.03630.05480.05290.0846
0.06800.05431.43740.07340.05350.05640.1000-0.3000-0.1000
0.10000.30000.700010.0000
3.00001.50002.5000Volatility forecast = 1.4868

29
TR4/2000

NP3533

5

ソフトウェアのテスト

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

初期パラメータ推定値は真値の半分になるようにします

GARCH(1,2) 数列のみモデル化されます

時間依存回帰ベクトル

定数要素、 1

線形ランプ要素、



正弦波要素,にはすべてのテストにおいて以下の形式をもつ要素があります：()


GARCH(p,q) 数列のモデリングの難しさは、p と q 両方に依存し、プロセス中にどのくらいボ
ラティリティメモリがあるかにも依存します。

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

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

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

30
TR4/2000

NP3533

6

モンテカルロの結果

AGARCH-type1
ここで使用される AGARCH(1,2)-type1 シミュレーションには以下のパラメータがあります：()

系列のショックは、4.10に設定された自由度のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、 (i 1…,3 )に続いて

、 、 、

、

次にdf （もしあれば）という降順に出力されます。

ガウス分布
推定値
0.22670.08730.14830.6944-0.2408-1.49822.5562-3.0455

推定標準誤差0.17120.05810.09030.10460.30240.25860.85950.6588推定値の標準誤差0.11470.06070.08160.07930.23980.24720.77280.5912

正しい値
0.20
0.10
0.15
0.70
-0.20
-1.50
2.50
-3.00

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

推定値
0.2022
0.09240.14680.7050-0.2109-1.50722.5011-3.0001

推定標準誤差0.06690.04110.05720.05170.15980.09620.16190.1695推定値の標準誤差
0.06390.0396
0.05260.04730.13800.09710.15710.1651

正しい値
0.20
0.10
0.15
0.70
-0.20
-1.50
2.50
-3.00

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

31
TR4/2000

NP3533

推定値
0.20380.09590.1526
0.6986-0.2076-1.5004
2.4993-2.9998

推定標準誤差0.04670.02850.03950.03570.09910.07080.0576
0.0952

推定値の標準誤差
0.04420.02830.03780.0332
0.09440.06610.05590.0902

正しい値
0.20
0.10
0.15
0.70
-0.20
-1.50
2.50
-3.00

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

スチューデントt-分布
推定値0.1983
0.09370.15280.6970-0.23214.3246-1.5000
2.4907-2.9950

推定標準誤差0.07060.0496
0.06840.0586
0.19540.66450.07650.12700.1331

推定値の標準誤差
0.0676
0.05920.07650.05820.18500.61890.08120.12740.1354

正しい値
0.20
0.10
0.15
0.70
-0.20
4.10
-1.50
2.00
-3.00

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

推定値
0.2011
0.09750.15230.6954-0.22034.2278-1.50642.4974-2.9987

推定標準誤差0.04630.03690.04850.03790.13380.4241
0.05240.04510.0743

推定値の標準誤差
0.0456
0.03370.04720.03940.12130.40690.05170.04390.0709

正しい値
0.20
0.10
0.15
0.70
-0.20
4.10
-1.50
2.00
-3.00

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

32
TR4/2000

NP3533

AGARCH-type2
ここで使用される AGARCH(1,2)-type2 シミュレーションには以下のパラメータがあります：()

系列のショックは、4.10に設定された自由度のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
The GARCH モデルのパラメータは、 (i 1…,3 )に続いて

、

、 、

、

次にdf （もしあれば）という降順に出力されます。

ガウス分布
推定値
0.12380.10200.14590.6778-0.1097-1.52452.4532-2.9564

推定標準誤差0.0651
0.06360.0847
0.09140.12980.19000.60220.4586

推定値の標準誤差
0.05720.06270.08560.08450.11160.17630.55500.4268

正しい値
0.10
0.10
0.15
0.70
-0.10
-1.50
2.50
-3.00

推定値の標準誤差
0.03150.0402
0.05380.0487
0.06130.06750.10890.1142

正しい値
0.10
0.10
0.15
0.70
-0.10
-1.50
2.50
-3.00

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

推定値
0.10810.10150.14820.6918-0.0993-1.50292.4971-2.9935

推定標準誤差0.03240.0436
0.05870.0441
0.06690.06780.10830.1116

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

33
TR4/2000

NP3533

推定値
0.10480.10040.14910.6960-0.1017-1.5006
2.4991-2.9996

推定標準誤差0.02190.03100.04020.03290.04610.04730.0356
0.0571

推定値の標準誤差
0.02150.02850.03790.03330.04250.04590.0391
0.0628

正しい値
0.10
0.10
0.15
0.70
-0.10
-1.50
2.50
-3.00

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

スチューデントt-分布
推定値0.1103
0.09840.15470.6856-0.1060
4.2912-1.4998
2.5087-3.0118

推定標準誤差0.03970.05480.0746
0.0666
0.0942
0.66390.04600.08530.0860推定値の標準誤差
0.03440.05680.07520.05980.08480.60370.05400.0882
0.0931

正しい値
0.10
0.10
0.15
0.70
-0.10
4.10
-1.50
2.00
-3.00

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

推定値0.10350.10050.15280.6924-0.0998
4.1854-1.49712.4986-2.9977

推定標準誤差0.0236
0.03600.04880.04170.05990.43450.03240.02840.0441

推定値の標準誤差
0.02210.03470.04840.04040.05690.39880.03600.03060.0494

正しい値
0.10
0.10
0.15
0.70
-0.10
4.10
-1.50
2.00
-3.00

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

34
TR4/2000

NP3533GJR-GARCH
ここで使用される GJR-GARCH(1,2) シミュレーションには以下のパラメータがあります：()

系列のショックは、4.10に設定された自由度数のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
The GARCH モデルのパラメータは、 (i 1…,3 )に続いて

、

、 、

、

次にdf （もしあれば）という降順に出力されます。

ガウス分布
推定値
0.11030.13670.20030.38050.1006-1.4983
2.4982-3.0025

推定標準誤差0.04100.07660.11850.14350.07660.10100.32880.2560推定値の標準誤差
0.03640.07740.10450.12700.07590.09920.31190.2397

正しい値
0.10
0.15
0.20
0.40
0.10
-1.50
2.50
-3.00

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

推定値
0.10480.14510.19970.3936
0.0958-1.5038
2.4975-2.9970

推定標準誤差0.02690.05370.0752
0.09270.05110.03790.06450.0677推定値の標準誤差
0.02210.0502
0.06690.07840.04780.03870.06270.0659

正しい値
0.10
0.15
0.20
0.40
0.10
-1.50
2.50
-3.0

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

35
TR4/2000

NP3533

推定値
0.10180.14700.20010.39720.1007-1.5016
2.4998-3.0000

推定標準誤差0.01720.03870.04990.06020.03850.0261
0.02290.0363

推定値の標準誤差
0.01520.03580.0476
0.05480.03400.02640.02240.0361

正しい値
0.10
0.15
0.20
0.40
0.10
-1.50
2.50
-3.00

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

スチューデントt-分布
推定値0.10270.13920.20470.39040.11584.2716-1.50132.4977-2.9989

推定標準誤差0.0266
0.0691
0.0932
0.0976
0.06550.66030.03280.0526
0.0549推定値の標準誤差
0.02850.08130.10280.11310.10670.63950.0336
0.05640.0599

正しい値
0.10
0.15
0.20
0.40
0.10
4.10
-1.50
2.00
-3.00

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

推定値
0.10080.14690.2016
0.39400.10674.2076-1.4982
2.4994-3.0003

推定標準誤差0.01780.05010.0716
0.06830.04890.43600.02210.01870.0299推定値の標準誤差
0.01700.0531
0.06850.06630.06060.40400.02310.01880.0314

正しい値
0.10
0.15
0.20
0.60
0.10
4.10
-1.50
2.00
-3.00

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

36
TR4/2000

NP3533

EGARCH
ここで使用される EGARCH(1,2) シミュレーションには以下のパラメータがあります：,()

系列のショックは、4.10に設定された自由度数のガウス分布かスチューデント t-分布のどちらかから得ています。
The GARCH モデルのパラメータは、 (i 1…,3 )に続いて

、

、

、

、

次にdf （もしあれば）という降順に出力されます。

ガウス分布
推定値
0.2786-0.1965-0.24810.06990.13460.30451.51812.56222.9430

推定標準誤差0.10160.07580.09580.11860.14450.22130.18550.64130.4929推定値の標準誤差
0.09740.07380.09090.12280.12850.18410.19280.6017
0.4626

正しい値
0.30
-0.20
-0.25
0.10
0.15
0.30
1.50
2.50
3.00

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

推定値
0.2958-0.1954-0.2533
0.08480.14320.2914
1.49702.49693.0004

推定標準誤差0.07090.04710.05950.07680.08500.14370.07470.12450.1301

推定値の標準誤差
0.06650.04620.05710.07890.08120.12470.07560.12310.1309

正しい値
0.30
-0.20
-0.25
0.15
0.10
0.30
1.50
2.50
3.00

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

37
TR4/2000

NP3533

推定値
0.2982-0.1980-0.2540
0.09320.14910.29671.50172.49992.9984

推定標準誤差0.04520.03170.03810.05080.05910.08950.05280.04420.0748推定値の標準誤差
0.04710.03260.0401
0.05330.05450.08700.05150.04310.0702

正しい値
0.30
-0.20
-0.25
0.15
0.10
0.30
1.50
2.50
3.00

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

スチューデントt-分布
推定値0.2811-0.2071-0.2452
0.08650.14640.32234.31191.49832.48713.0094

推定標準誤差0.10220.06280.0786
0.09290.09320.17380.68190.06680.11420.1195推定値の標準誤差
0.09240.06390.07940.09030.09240.1526
0.41590.07380.12100.1267

正しい値
0.30
-0.20
-0.25
0.10
0.15
0.30
4.10
1.50
2.50
3.00

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

推定値
0.2859-0.2016-0.2483
0.09260.15120.31544.22811.49422.49763.0001

推定標準誤差0.07020.0426
0.05940.0652
0.06330.11740.43480.04570.04130.0681

推定値の標準誤差
0.08110.04930.0706
0.06570.08010.12830.38540.05390.04510.0734

正しい値
0.30
-0.20
-0.25
0.10
0.15
0.30
4.10
1.50
2.50
3.00

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

38
TR4/2000

NP3533

結論
表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 関数の呼び出しのほうがソフトウェア開発者にとって簡単に行うことができると言えます。時系列パッケージは対話型インターフェース39
TR4/2000

NP3533

の利点がありますが、大規模システムの構成要素として時系列パッケージをプログラムで呼び出すのは容易ではありません（あるいは計算効率が良くはありません）。現在の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 ofUnited 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 theVolatility 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

40
TR4/2000

NP3533