統計ソフト
R
によるボラティリティ変動モデル入門
∗
松本 惇
†1
はじめに
金融の世界において、ボラティリティとは、間単に言えば、ある金融商品の分散又は標準偏差の事であ る．このボラティリティは、投資リスクを表す一つの指標であり、オプションの世界においては、その価 格を決定する要因の一つとなっている. このボラティリティが、時間を通じて変化するのであれば、その 変化の特徴を理解し、予想する事は、そのようなリスクを管理する上で非常に重要となろう. 本稿では、 統計ソフト R を用いる事で、様々なボラティリティ変動モデルを推定し、将来のボラティリティ変動を 予測する方法を述べる. 統計ソフト R（以下、R）では、時系列分析を行うためのパッケージとして、tseries と fSeries の 二種類がポピュラーであり、このいずれを用いても基本的なボラティリティ変動モデルを推定したり、予 測したりする事ができる1．まず、tseries と fSeries をインストールする. R を起動して、ツールバー に「パッケージ」とあるので、選択して、次に「パッケージのインストール」を選択すれば良い. ダウ ンロード先が、多数表示されるかもしれないが、とりあえず Japan(Tokyo) あたりを選び、パッケージ の名前を選択する. インストールが完了すれば、R の画面上で library(tseries)、library(fSeries) と入力して、使 えるようにしておかなければならない.2
ボラティリティ変動モデルの推定
ある金融商品の t 時点における収益率 rtは、t − 1 時点において利用可能な情報集合 Ωt−1により予測 される変動 E[rt|Ωt−1] と、予測不可能なショック ²tに分割する事ができる. rt= E[rt|Ωt−1] + ²t.ここで、²tは次のように表されるとする2. 但し、iid は独立同一分布 (independently, identically
dis-tributed) を意味する．
²t= σtzt, σt> 0, zt∼ iid N (0, 1). ∗推定には、統計ソフト R Ver.2.4.1 を使用した.
1fSeries の方が色々と便利だが、信頼性の面で tseries に劣る. また、fSeries で GARCH モデルを推定する際には、
garchFit という関数を用いるが、この関数が正常に動作しない状況は少なくない．
2渡部 [2000] では、z
tは独立同一分布に従うとされ、正規分布に従うと言われていない．しかし、後述のように最尤法で推定 を行う際には、ztが正規分布に従うと仮定するため、本稿では、この時点で ztの正規性を仮定した．
我々の関心は、σtであり、この σt又は、σ2t をボラティリティと呼ぶ． 本稿では、ボラティリティ変動モデルの推定に際し、日経平均 225 の日次データを用いる．原系列をそ のまま推定に用いる事はできないため、その処理も含めて、以下、推定方法を述べる． 使用するデータを nikkei と表記する．これは、大屋先生の HP にある私のページからダウンロードで きる．実際のファイルには、nikkei4.csv と名前がついているが、ダウンロードして、H ドライブに保 存する．その後、R でこのファイルを読み込むためには、 nikkei4<-read.csv("h:=Y=Ynikkei4.csv",header=T) (1) とすれば良い．nikkei の一行目は、データ系列の名前であり、その名前も同時に読み込むため、header=T としている．この nikkei4.csv の 4 列目に、終値が含まれているので、ここでは、その終値を用いる． 4 列目だけを取り出すためには、 nikkei4<-data.frame(nikkei4) (2) としてから、終値の系列の名前を nikkei として、 nikkei<-nikkei4[,4] (3) とすれば良い．我々の関心は、実際の日経平均の値ではなく、その日次変化率である．変化率は、対数 差分で近似する事ができるため、 d.nikkei <-diff(log(nikkei)) (4) と変換する事で、新しい系列 d.nikkei を生成する． 以下、この d.nikkei を使って、説明を続ける．
2.1
ARCH(q) モデルの推定
ARCH(:Auto-regressive Conditional Heteroskedasticity) は、ボラティリティσ2
t を過去の収益率の予 測不可能なショックの二乗 ²2t−1, ²2 t−2, · · · の線形関数をして次のように定式化される． ARCH(q) : σ2t = ω + q X j=1 αj²2t−j, ω > 0, ∀j αj ≥ 0.
tseries において、ARCH モデルを推定するためには、garch コマンドを使う．例えば、ARCH(2) モ デルを推定したいのであれば、推定モデルの名前を arch2 として、次のように入力する
arch2<-garch(d.nikkei, order=c(0,2)) (5) 次に、summary(arch2) と入力すれば、次のような出力がなされる3.
3***** ESTIMATION WITH ANALYTICAL GRADIENT *****と表示されれば、パラメータの収束に関しての問題は無い．もし、
Call:
garch(x = d.nikkei, order = c(0, 2)) Model:
GARCH(0,2) Residuals:
Min 1Q Median 3Q Max -7.55848 -0.54383 0.01040 0.56867 6.09377 Coefficient(s):
Estimate Std. Error t value Pr(>|t|) a0 1.366e-04 2.783e-06 49.077 < 2e-16 *** a1 1.157e-01 1.452e-02 7.971 1.55e-15 *** a2 8.335e-02 1.487e-02 5.606 2.07e-08 ***
---Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Diagnostic Tests:
Jarque Bera Test data: Residuals
X-squared = 1611.481, df = 2, p-value < 2.2e-16 Box-Ljung test data: Squared.Residuals X-squared = 1.1212, df = 1, p-value = 0.2897 推定された ARCH(2) モデルは、次のようになる． σ2t = 1.366 × 10−4+ 0.1157²2t−1+ 0.08335²2t−2, Std. Error: (2.783 × 10−6) (0.01452) (0.01487).
2.2
GARCH(p, q) モデルの推定
GARCH(:Generalized Auto-regressive Conditional Heteroskedasticity) は、ボラティリティσ2
tを過去 の収益率の予測不可能なショックの二乗 ²2t−1, ²2 t−2, · · · と過去のボラティリティσ2t−1, σt−22 , · · · の線形関 数をして次のように定式化される． GARCH(p, q) : σ2 t = ω + p X i=1 βiσt−12 + q X j=1 αj²2t−j, ω > 0, ∀i βi ≥ 0, ∀j αj ≥ 0.
GARCH モデルを推定する際にも、同じく garch コマンドを用いる．例えば、GARCH(1,1) を推定する のであれば、推定モデルの名前を garch11 として、
garch11<-garch(d.nikkei, order=c(1,1)) (6) と入力し、続いて、summary(garch11) とすれば、次のような出力がなされる．
Call:
garch(x = d.nikkei, order = c(1, 1)) Model:
GARCH(1,1) Residuals:
Min 1Q Median 3Q Max -8.98948 -0.58079 0.01014 0.59228 5.88422 Coefficient(s):
Estimate Std. Error t value Pr(>|t|) a0 3.886e-06 4.801e-07 8.094 6.66e-16 ***
a1 6.483e-02 4.907e-03 13.212 < 2e-16 *** b1 9.134e-01 6.331e-03 144.286 < 2e-16 ***
---Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Diagnostic Tests:
Jarque Bera Test data: Residuals
X-squared = 2179.503, df = 2, p-value < 2.2e-16 Box-Ljung test data: Squared.Residuals X-squared = 1.0566, df = 1, p-value = 0.304 この場合、次の GARCH(1,1) モデルが推定された事になる． σ2 t = 3.886 × 10−6+ 0.06483²2t−1+ 0.9134σt−12 , Std. Error: (4.801 × 10−7) (0.0049) (0.006331).
2.3
ARCH,GARCH モデルの診断
ここでは、推定した GARCH(1,1) モデルの推定結果について、診断を行う．summary(garch11) と入 力すれば、GARCH(1,1) モデルの推定結果が出力されるが、その中に、 Diagnostic Tests:Jarque Bera Test data: Residuals
X-squared = 2179.503, df = 2, p-value < 2.2e-16 Box-Ljung test
data: Squared.Residuals
X-squared = 1.0566, df = 1, p-value = 0.304
と書かれている．tseries では、GARCH モデルを推定した際、この二つの検定から、モデルの当ては まりについての診断を行う事ができる．
Jarque Bera Test: この検定の帰無仮説は、「残差系列が正規分布に従う」である．残差系列が正規 性を持っている事が望ましいが、このモデルにおいては、その p 値が非常に小さいため、帰無仮説を棄 却する事になる． Ljung-Box Test: この検定の帰無仮説は、「残差の二乗に 1 次の系列相関がある」である4．そのよう な系列相関は無い方が望ましいが、このモデルでは、p 値が大きく、帰無仮説を採択する． また、視覚的に診断を行う事もできる．次のプログラム5を入力すれば、下の図のような結果が得られる． n<-length(d.nikkei) nparam<-length(coef(garch11))-1 L<-30 garch.box<-rep(0,L) garch.boxsq<-rep(0,L) for(i in 1:L){ box<-Box.test(garch11$resid,lag=i,type="Ljung-Box") boxsq<-Box.test(garch11$resid^2,lag=i,type="Ljung-Box") garch.box[i]<-box$p.value 4残差ではなく、残差の二乗を考えるのは、GARCH が ²tをモデル化したものではなく、その分散をモデル化しているからで ある．
5このプログラムは、大屋先生の HP にある私のページからダウンロードできる．なお、acf や pacf の所で、plot=FALSE と
garch.boxsq[i]<-boxsq$p.value } par(mfrow=c(3,4),mar=c(4,4,4,0.5),las=1) plot(garch11$resid,type="l",main="Residual") plot(garch.box,type="b",ylim=c(0,1),main="Ljung-Box p-value") abline(h=0.05,lty=2,col="red") abline(v=2.5) acf(garch11$resid[nparam:n],main="ACF") pacf(garch11$resid[nparam:n],main="PACF") cpgram(garch11$resid,main="Cumulative Periodgram") qqnorm(garch11$resid,main="Normal QQ plot") qqline(garch11$resid) plot(garch11$resid,type="n",bty="n",xlab="",ylab="",xaxt="n",yaxt="n") plot(garch11$resid^2,type="l",main="Residual Square") plot(garch.boxsq,type="b",ylim=c(0,1),main="Ljung-Box p-value") abline(h=0.05,lty=2,col="red") abline(v=2.5) acf(garch11$resid[nparam:n]^2,main="ACF") pacf(garch11$resid[nparam:n]^2,main="PACF") 今、３行４列に渡り、図がプロットされているが、各々、次の事を示している． • 1 行 1 列目：残差系列の時系列変動を表している． • 1 行 2 列目：Ljung-Box 検定における p 値を示している． 即ち、残差系列に (30 次までの) 系列相関が無いかどうかを検定した場合の p 値を示している．よっ て、この p 値は大きいほうが良い． • 1 行 3 列目：残差の ACF（自己相関関数）を（35 次まで) 示している． 破線は、5%の信頼度で 0 である範囲を示しているので、この破線内に含まれている方が良い. • 1 行 4 列目：残差の PACF(偏自己相関関数）を（35 次まで）示している. 解釈については、上と同じ． • 2 行 1 列目：累積ピリオドグラムを示している． 破線が平行に二本並んでいるが、その間に実線が収まっていれば、残差は正規性を持つ．しかし、 ここでの図では、大きく外れているため、正規性を持つとは言えない． • 2 行 2 列目：正規 QQ プロットを示している． 中央に走る実線上に太い線が位置していれば、残差は正規性を持つ． • 2 行 4 列目：残差二乗系列の時系列変動を表している． • 3 行 1 列目：Ljung-Box 検定における p 値を示している． ここでは、残差二乗系列に（30 次までの）系列相関が無いかを検定した場合の p 値を示している． • 3 行 2 列目：残差二乗系列の ACF を（35 次まで) 示している． • 3 行 3 列目：残差二乗系列の PACF を（35 次まで）示している．
3
応用的な内容
3.1
EGARCH(p, q) モデルと GJR(p, q) モデルの推定
前述の ARCH モデルや GARCH モデルを、例えば株式収益率のボラティリティの分析に応用する事 を考える．その際に我々が直面する問題は、ボラティリティ変動の非対称性を ARCH モデルや GARCH モデルが捉える事が出来ないという点である．そのため、ボラティリティ変動の非対称性を捉える事が 出来るモデルとして、ここでは、EGARCH モデルと GJR モデルを用いる． GJR(p, q) モデル: GJR(p, q) モデルは、²t−1が負であれば 1、それ以外には 0 となるようなダミー変数 D−t−1を用いて、ボラティリティ変動を次のように定式化する． σt2= ω + p X i=1 βiσ2t−i+ q X j=1 (αj²2t−j+ γjD−t−j²2t−j), ω > 0, ∀i, j βi, αj, γj≥ 0, . EGARCH(p, q) モデル: これまでのモデルでは、ボラティリティσ2 t の説明変数として、過去のボラ ティリティや、ショックの二乗 ²2t−jなどを考えていた．しかし、EGARCH モデルでは、²t−jではなく、 それをボラティリティの平方根である σt−jで割り、基準化した zt−j := ²t−j/σt−jを説明変数として用 い、次のようにボラティリティ変動を定式化する． ln(σ2 t) = ω + p X i=1 βiln(σt−i2 ) + q X j=1 αj © θzt−j+ γ(|zt−j| − E|zt−j|) ª . これらのモデルを R を用いて推定する際には、fSeries に含まれている garchOxInterface という関 数を使う．しかし、単に fSeries をインストールしただけでは、この関数を使う事は出来ない．以下、 garchOxInterface を作動させるためのステップを説明する． Step1: Ox をインストールする．Ox は、例えば、以下のサイトからダウンロードできる．インストー ル先は (経済学研究科の PC ルームなら)、H:=YOxMetrics4=YOx とする. その後、OxMetrics フォル ダにある Ox フォルダを、コピーして、H:=YOx にコピーする6.
東京大学大森裕浩先生のサイト : http://www.e.u-tokyo.ac.jp/∼omori/Ox/index.htm
Step2: Ox をインストールした後に、Ox のパッケージである [email protected] をダウンロードする．こ れは、以下のサイトからダウンロードできるが、その際、自分のメールアドレスを登録しなければ ならない．登録したアドレスに、ダウンロードサイトのアドレスが送られてくる. ダウンロードし た後は、H:=YOx=Ypackages で解凍する. このフォルダが無い場合には作成する.
http://www.core.ucl.ac.be/∼laurent/G@RCH/site/index.html
Step3 : GarchOxModeling.ox というファイルを以下のサイトからダウンロードする. このサイトのペー ジ中程に、Modified GarchOxModelling file と書かれている. ダウンロードして、H:=YOx=Ylib にコピーする.
http://faculty.chicagogsb.edu/ruey.tsay/teaching/bs41202/sp2007/
Step4 : Step3 と同じサイトから、A corrected version of the garchOxFit function. をダウン ロードする．先程ダウンロードした際とであれば、、Modified GarchOxModelling file のすぐ 下に掲載されている．これは、txt 形式で保存されているので、名前をつけて保存する．その際、 名前が
garchOxFit R.txt
となっている事を確認する．その後で、R の作業フォルダにコピーする7．次に、このファイルを 修正する．ファイルを開いて、一番下から 16 行目に、
command = "C:=Y=YOx=Y=Ybin =Y=Yoxl.exe C:=Y=YOx=Y=Ylib=Y=YGarchOxModelling.ox"
と書かれているので、C:=Y の部分を、今なら、H:=Y と書き換える．Ox を H ドライブ以外にインス トールしたなら、そのフォルダの名前にする．最後に、R 上で source("garchOxFit R.txt") と入力して、エラーが出なければ、推定のための準備が完了する． モデルの推定: garchOxInterface を使うと、単にボラティリティ変動モデルを当てはめるだけでなく、 使用するデータ系列に ARMA モデルを当てはめ、その後にボラティリティ変動モデルを当てはめる． よって、平均が ARMA(1,0)、即ち、AR(1) に従い、ボラティリティが GJR(1,1) に従うなら、そのモデ ルの名前を m1 として、
m1<-garchOxFit(formula.mean=∼arma(1,0), formula.var=∼gjr(1,1),series=dif) (7)
とすれば良い．その結果、次のような結果が得られる8. ********************
** SPECIFICATIONS ** ******************** Dependent variable : X
Mean Equation : ARMA (1, 0) model. No regressor in the mean
Variance Equation : GJR(1, 1) model. No regressor in the variance
The distribution is a Gauss distribution. Strong convergence using numerical derivatives Log-likelihood = 11698.9
Maximum Likelihood Estimation (Std.Errors based on Second derivatives) Coefficient Std.Error t-value t-prob
Cst(M) 0.000101 0.00021415 0.4725 0.6366 AR(1) 0.150090 0.017042 8.807 0.0000 Cst(V) 0.040843 0.0080254 5.089 0.0000 ARCH(Alpha1) 0.068120 0.0074708 9.118 0.0000 GARCH(Beta1) 0.908588 0.0099891 90.96 0.0000 GHR(Gamma1) 0.090895 0.0336651 2.700 0.0071 No. Observations : 3920 No. Parameters : 6 Mean (Y) : -0.00010 Variance (Y) : 0.00017 Skewness (Y) : -0.16499 Kurtosis (Y) : 5.99703 Log Likelihood : 11698.930 Alpha[1]+Beta[1]: 0.97651 Warning : To avoid numerical problems, the estimated parameter Cst(V), and its std.Error have been multiplied by 10^4.
The sample mean of squared residuals was used to start recursion. The positivity constraint for the GARCH (1,1) is observed.
This constraint is alpha[L]/[1 - beta(L)] >= 0.
7作業フォルダを確認するには、R を起動して、「ファイル」→「ディレクトリの変更」と進む．これで、作業フォルダのアドレ スが表示されるので、その場所に garchOxFit R.txt をコピーすれば良い．また、R のインストール先は、デフォルトのままで ある事．デフォルトでは、R は C ドライブの Program Files 内にインストールされるが、これを変更してインストールすると、 garchOxFit はエラーを起こす． 8推定には、GARCH モデルと比べて、非常に時間がかかる．CPU が 1.5GHz、メモリーが 512MB 程度では、数分かかる場 合があり、ノート型 PC では負荷のため、フリーズする場合もある．
The unconditional variance is 0.000173863
The conditions are alpha[0] > 0, alpha[L] + beta[L] < 1 and alpha[i] + beta[i] >= 0. => See Doornik & Ooms (2001) for more details.
The condition for existence of the fourth moment of the GARCH is observed. The constraint equals 0.962849 and should be < 1.
=> See Ling & McAleer (2001) for details. Estimated Parameters Vector :
0.000101; 0.150090; 0.000004; 0.068120; 0.908588 ***************
** FORECASTS ** ***************
Number of Forecasts: 15 Horizon Mean Variance
1 0.002854 0.0001076 2 0.0005144 0.000105 3 0.0001632 0.0001026 4 0.0001105 0.0001002 5 0.0001026 9.78e-005 6 0.0001014 9.55e-005 7 0.0001012 9.326e-005 8 0.0001012 9.107e-005 9 0.0001012 8.893e-005 10 0.0001012 8.684e-005 11 0.0001012 8.48e-005 12 0.0001012 8.281e-005 13 0.0001012 8.087e-005 14 0.0001012 7.897e-005 15 0.0001012 7.711e-005
---Forecasts errors measures cannot be computed because there are not enough out-of-sample observations).
Elapsed Time : 4.031 seconds (or 0.0671833 minutes). ***********
** TESTS ** ***********
Statistic t-Test P-Value Skewness -0.38815 9.9251 3.2391e-023 Excess Kurtosis 3.5788 45.767 0.00000 Jarque-Bera 2190.3 .NaN 0.00000
---Information Criterium (to be minimized)
Akaike -5.966291 Shibata -5.966294 Schwarz -5.958288 Hannan-Quinn -5.963451
---Q-Statistics on Standardized Residuals
--> P-values adjusted by 1 degree(s) of freedom Q( 10) = 8.69102 [0.4662708]
Q( 15) = 16.6551 [0.2750242] Q( 20) = 25.9216 [0.1323881]
H0 : No serial correlation ==> Accept H0 when prob. is High [Q < Chisq(lag)]
---Q-Statistics on Squared Standardized Residuals --> P-values adjusted by 2 degree(s) of freedom Q( 10) = 6.63530 [0.5764419]
Q( 15) = 8.15385 [0.8334345] Q( 20) = 11.0363 [0.8928108]
H0 : No serial correlation ==> Accept H0 when prob. is High [Q < Chisq(lag)]
---ARCH 1-2 test: F(2,3913)= 0.86756 [0.4201] ARCH 1-5 test: F(5,3907)= 0.74896 [0.5868] ARCH 1-10 test: F(10,3897)= 0.57124 [0.8387]
---Diagnostic test based on the news impact curve (EGARCH vs. GARCH) Test P-value
Sign Bias t-Test 0.73704 0.46110 Negative Size Bias t-Test 0.14646 0.88356 Positive Size Bias t-Test 1.25593 0.20914 Joint Test for the Three Effects 1.64696 0.64879
---Joint Statistic of the Nyblom test of stability: 27.826 Individual Nyblom Statistics:
Cst(M) 0.12959 AR(1) 0.46767 Cst(V) 20.91820 ARCH(Alpha1) 25.98398 GARCH(Beta1) 22.27730 GJR(Gamma1) 20.11231
Rem: Asymptotic 1% critical value for individual statistics = 0.75. Asymptotic 5% critical value for individual statistics = 0.47.
---Adjusted Pearson Chi-square Goodness-of-fit test
# Cells(g) Statistic P-Value(g-1) P-Value(g-k-1) 40 128.5306 0.000000 0.000000 50 145.9184 0.000000 0.000000 60 161.8061 0.000000 0.000000 Rem.: k = 5 = # estimated parameters
---この推定結果から、次のモデルが推定された事が分かる． rt= 0.000101 + 0.15009rt−1+ ²t, Std. Error: (0.000214) (0.017042), σ2 t = 0.040843 + 0.908588σt−12 + (0.06812 + 0.090895Dt−1− )²2t−1, Std. Error: (0.00802) (0.00998) (0.0076) (0.0336). また、平均が AR(1) に従い、ボラティリティが EGARCH(1,1) に従うなら、そのモデルを m2 として、 m2<-garchOxFit(formula.mean=∼arma(1,0), formula.var=∼egarch(1,1),series=dif) とすれば良い．
3.2
モデルにおける次数の決定
実際にモデルを推定する際には、モデルの次数を決める必要がある．その際には、AIC を用いて決定す る事が便利である．しかし、逐次モデルを推定していき、AIC を求めるのは面倒であり、しかも、garch コマンドでは、AIC を自動的に出力しない．よって、ここでは、次のような手順で求める．まず、AIC は次のように定義される． AIC := −2 log bL + 2n. 但し、log bL は、推定されたパラメータの下で計算された対数尤度であり、n は推定されたパラメータ数 である．GARCH(p, q) の場合には、n = p + q + 1 となる．この時、次のようなプログラムを考える9． for(p in 1:3){ for(q in 1:3){ tmp<-garch(dif,order=c(p,q)) likelihood<-tmp$n.likeli 9このプログラムは、大屋先生の HP にある私のページからダウンロードできる.nparam<-length(coef(tmp)) aic<- -2*likelihood+2*nparam
cat("order1",p,"order2",q,"AIC",aic, "\n")}}
これにより、9 種類のモデルを自動的に推定し、AIC を出力する．
4
参考資料
Doornik, J.A. (2002), Object-Oriented Matrix Programming Using Ox, 3rd ed. London: Timber-lake Consultants Press and Oxford: www.doornik.com.
渡部敏明 (2004), ボラティリティ変動モデル, 朝倉書店 また、以下のサイトから、tseries と fSeries のマニュアルをダウンロードできる． http://cran.r-project.org/src/contrib/Descriptions/tseries.html http://cran.r-project.org/src/contrib/Descriptions/fSeries.html 加えて、Ox のインストールや、garchOxFit の使用準備については、次のサイトが参考になる． http://faculty.chicagogsb.edu/ruey.tsay/teaching/bs41202/sp2007/
