関連情報

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  参考文献

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

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

1  概要

時間依存性のある分散の数列のモデリングは金融工学の多くの分野で重要です。このテクニカルレポートでは、NAG社で数値計算ライブラリ向けに開発された、一変量の一般化自己回帰条件つき分散不均一(GARCH)ルーチンについて記述しています。ここではFortran 77 ライブラリ [1] の次期バージョン用に開発されたGARCH ソフトウェアに扱いますが、GARCH ソフトウェアのいくつかは既に現バージョンのNAG C ライブラリ [2] に含まれています。ここで記述しているFortran 77 ソフトウェアは任意の値 pq をもつ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] は以下の形式で表されます:

figure1
figure2

このプロセスは、q + 1 個の係数 αii = 0,…,qp 個の係数 βii = 1,…,p、平均 b0k 個の線形回帰係数bii = 1,…,k、内生変数 yt /外生変数 xt 、ショックεt 、条件付き分散 ht 、そして時間 t までの一連の情報 Ψt で記述されます。
p = 0 の場合に GARCH (p,q) モデルはARCH (q) モデルとも呼ばれます。

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

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

AGARCH (p,q)-type1

figure3

AGARCH (p,q)-type2

figure4

GJR-GARCH (p,q), または Glosten-Jagannathan-Runkle GARCH [7]
figure5
ここでは εt > 0 ならば St = 1 また εt ≥ 0 ならばSt = 0

EGARCH (p,q), または指数 GARCH
figure6
ここでは、figure7E[| Zt-i |] は | Zt-i |の期待値を表します。

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

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

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

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

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

上記の全ての GARCH プロセスは、パラメータベクトルθで独自に記述されます。 この場合、 θ = (b0,bT,ωT), ωT = (α0, α1, α2,…,αq,β 1,β2,…,βp,γ) かつ bT = (b1,…,bk)です。

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

figure9

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

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

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

figure10

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 数列の条件つき分散, hi (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 配列.
呼び出し時、Xi 番目の行は時間依存の外生変数ベクトルを含みます。 xi (i = 1,…,T) 

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

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

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

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

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

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

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 配列.
復帰時、推定条件つき分散, hi (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.
呼び出し時、モデル化された数列からの、配列 HTET の項数

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)プロセスの過去の条件つき分散の数列, hi (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モデル、yt = b0 + εt
  • スチューデントt-分布からのショックと観測値をもつAGARCH (1,2)-type2モデル、 figure11

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


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モデル、figure11
  • スチューデントt-分布からのショックと観測値をもつAGARCH (1,1)-type2モデル, figure11

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


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)モデル、figure11
  • スチューデントt-分布からのショックと観測値をもつGJR-GARCH (1,1)モデル, figure11

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


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)モデル、figure11
  • スチューデントt-分布からのショックと観測値をもつEGARCH (1,2)モデル, figure11

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


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) 数列のみモデル化されます

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

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

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

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

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

6  モンテカルロの結果

AGARCH-type1

ここで使用される AGARCH (1,2)-type1 シミュレーションには以下のパラメータがあります:
k = 3
α0 = 0.2, α1 = 0.1, α2 = 0.15, β1 = 0.7
γ = -0.2
b0 = 0, b1 = -1.5, b2 = 2.5, b3 = -3.0
figure14
系列のショックは、4.10に設定された自由度のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、bi(i = 1,…,3 )に続いて α0α1, α2β1γ 次にdf (もしあれば)という降順に出力されます。

ガウス分布
推定値推定標準誤差推定値の標準誤差正しい値
0.22670.17120.11470.20
0.08730.05810.06070.10
0.14830.09030.08160.15
0.69440.10460.07930.70
-0.24080.30240.2398-0.20
-1.49820.25860.2472-1.50
2.55620.85950.77282.50
-3.04550.65880.5912-3.00
表 1::400 個の観測値と300 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.20220.06690.06390.20
0.09240.04110.03960.10
0.14680.05720.05260.15
0.70500.05170.04730.70
-0.21090.15980.1380-0.20
-1.50720.09620.0971-1.50
2.50110.16190.15712.50
-3.00010.16950.1651-3.00
表 2:1000 個の観測値と 200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.20380.04670.04420.20
0.09590.02850.02830.10
0.15260.03950.03780.15
0.69860.03570.03320.70
-0.20760.09910.0944-0.20
-1.50040.07080.0661-1.50
2.49930.05760.05592.50
-2.99980.09520.0902-3.00
表 3:2000 個の観測値と 200 回のシミュレーション

スチューデントt-分布
推定値推定標準誤差推定値の標準誤差正しい値
0.19830.07060.06760.20
0.09370.04960.05920.10
0.15280.06840.07650.15
0.69700.05860.05820.70
-0.23210.19540.1850-0.20
4.32460.66450.61894.10
-1.50000.07650.0812-1.50
2.49070.12700.12742.00
-2.99500.13310.1354-3.00
表 4:1000個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.20110.04630.04560.20
0.09750.03690.03370.10
0.15230.04850.04720.15
0.69540.03790.03940.70
-0.22030.13380.1213-0.20
4.22780.42410.40694.10
-1.50640.05240.0517-1.50
2.49740.04510.04392.00
-2.99870.07430.0709-3.00
表 5:2000 個の観測値と 200 回のシミュレーション

AGARCH-type2

ここで使用される AGARCH (1,2)-type2 シミュレーションには以下のパラメータがあります:
k = 3
α0 = 0.1, α1 = 0.1, α2 = 0.15, β1 = 0.7
γ = -0.1
b0 = 0, b1 = -1.5, b2 = 2.5, b3 = -3.0
figure14
系列のショックは、4.10に設定された自由度のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、bi(i = 1,…,3 )に続いて α0α1, α2β1γ 次にdf (もしあれば)という降順に出力されます。

ガウス分布
推定値推定標準誤差推定値の標準誤差正しい値
0.12380.06510.05720.10
0.10200.06360.06270.10
0.14590.08470.08560.15
0.67780.09140.08450.70
-0.10970.12980.1116-0.10
-1.52450.19000.1763-1.50
2.45320.60220.55502.50
-2.95640.45860.4268-3.00
表 6:400 個の観測値と300 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.10810.03240.03150.10
0.10150.04360.04020.10
0.14820.05870.05380.15
0.69180.04410.04870.70
-0.09930.06690.0613-0.10
-1.50290.06780.0675-1.50
2.49710.10830.10892.50
-2.99350.11160.1142-3.00
表 7:1000 個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.10480.02190.02150.10
0.10040.03100.02850.10
0.14910.04020.03790.15
0.69600.03290.03330.70
-0.10170.04610.0425-0.10
-1.50060.04730.0459-1.50
2.49910.03560.03912.50
-2.99960.05710.0628-3.00
表 8:2000 個の観測値と200 回のシミュレーション

スチューデントt-分布
推定値推定標準誤差推定値の標準誤差正しい値
0.11030.03970.03440.10
0.09840.05480.05680.10
0.15470.07460.07520.15
0.68560.06660.05980.70
-0.10600.09420.0848-0.10
4.29120.66390.60374.10
-1.49980.04600.0540-1.50
2.50870.08530.08822.00
-3.01180.08600.0931-3.00
表 9:1000 個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.10350.02360.02210.10
0.10050.03600.03470.10
0.15280.04880.04840.15
0.69240.04170.04040.70
-0.09980.05990.0569-0.10
4.18540.43450.39884.10
-1.49710.03240.0360-1.50
2.49860.02840.03062.00
-2.99770.04410.0494-3.00
表 10:2000 個の観測値と200 回のシミュレーション

GJR-GARCH

ここで使用される GJR-GARCH (1,2) シミュレーションには以下のパラメータがあります:
k = 3
α0 = 0.1, α1 = 0.15, α2 = 0.2, β1 = 0.4
γ = 0.1
b0 = 0, b1 = -1.5, b2 = 2.5, b3 = -3.0
figure14
系列のショックは、4.10に設定された自由度数のガウス分布あるいはスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、bi(i = 1,…,3 )に続いて α0α1, α2β1γ 次にdf (もしあれば)という降順に出力されます。

ガウス分布
推定値推定標準誤差推定値の標準誤差正しい値
0.11030.04100.03640.10
0.13670.07660.07740.15
0.20030.11850.10450.20
0.38050.14350.12700.40
0.10060.07660.07590.10
-1.49830.10100.0992-1.50
2.49820.32880.31192.50
-3.00250.25600.2397-3.00
表 11:400 個の観測値と300 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.10480.02690.02210.10
0.14510.05370.05020.15
0.19970.07520.06690.20
0.39360.09270.07840.40
0.09580.05110.04780.10
-1.50380.03790.0387-1.50
2.49750.06450.06272.50
-2.99700.06770.0659-3.00
表 12:1000 個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.10180.01720.01520.10
0.14700.03870.03580.15
0.20010.04990.04760.20
0.39720.06020.05480.40
0.10070.03850.03400.10
-1.50160.02610.0264-1.50
2.49980.02290.02242.50
-3.00000.03630.0361-3.00
表 13:2000 個の観測値と 200 回のシミュレーション

スチューデントt-分布
推定値推定標準誤差推定値の標準誤差正しい値
0.10270.02660.02850.10
0.13920.06910.08130.15
0.20470.09320.10280.20
0.39040.09760.11310.40
0.11580.06550.10670.10
4.27160.66030.63954.10
-1.50130.03280.0336-1.50
2.49770.05260.05642.00
-2.99890.05490.0599-3.00
表 14:1000 個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.10080.01780.01700.10
0.14690.05010.05310.15
0.20160.07160.06850.20
0.39400.06830.06630.60
0.10670.04890.06060.10
4.20760.43600.40404.10
-1.49820.02210.0231-1.50
2.49940.01870.01882.00
-3.00030.02990.0314-3.00
表 15:2000 個の観測値と 200 回のシミュレーション

EGARCH

ここで使用される EGARCH (1,2) シミュレーションには以下のパラメータがあります:
k=3, α0 = 0.3, α1 = -0.2, α2 = -0.25,
Φ1 = 0.1, Φ2 = 0.15, β1 = 0.3
b0 = 0, b1 = 1.5, b2 = 2.5, b3 = 3.0
figure15
系列のショックは、4.10に設定された自由度数のガウス分布かスチューデント t-分布のどちらかから得ています。
GARCH モデルのパラメータは、bi(i = 1,…,3 )に続いて α0α1, α2β1γ 次にdf (もしあれば)という降順に出力されます。

ガウス分布
推定値推定標準誤差推定値の標準誤差正しい値
0.27860.10160.09740.30
-0.19650.07580.0738-0.20
-0.24810.09580.0909-0.25
0.06990.11860.12280.10
0.13460.14450.12850.15
0.30450.22130.18410.30
1.51810.18550.19281.50
2.56220.64130.60172.50
2.94300.49290.46263.00
表 16:400 個の観測値と300 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.29580.07090.06650.30
-0.19540.04710.0462-0.20
-0.25330.05950.0571-0.25
0.08480.07680.07890.15
0.14320.08500.08120.10
0.29140.14370.12470.30
1.49700.07470.07561.50
2.49690.12450.12312.50
3.00040.13010.13093.00
表 17:1000 個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.29820.04520.04710.30
-0.19800.03170.0326-0.20
-0.25400.03810.0401-0.25
0.09320.05080.05330.15
0.14910.05910.05450.10
0.29670.0895 0.08700.30
1.5017 0.05280.05151.50
2.49990.04420.0431 2.50
2.99840.07480.07023.00
表 18:2000 個の観測値と200 回のシミュレーション

スチューデントt-分布
推定値推定標準誤差推定値の標準誤差正しい値
0.28110.10220.09240.30
-0.20710.06280.0639-0.20
-0.24520.07860.0794-0.25
0.08650.09290.09030.10
0.14640.09320.09240.15
0.32230.17380.15260.30
4.31190.68190.41594.10
1.49830.06680.07381.50
2.48710.11420.12102.50
3.00940.11950.12673.00
表 19:1000 個の観測値と200 回のシミュレーション

推定値推定標準誤差推定値の標準誤差正しい値
0.28590.07020.08110.30
-0.20160.04260.0493-0.20
-0.24830.05940.0706-0.25
0.09260.06520.06570.10
0.15120.06330.08010.15
0.31540.11740.12830.30
4.22810.43480.38544.10
1.49420.04570.05391.50
2.49760.04130.04512.50
3.00010.06810.07343.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

Results matter. Trust NAG.

Privacy Policy | Trademarks