本文／ＹＡ６２９６Ｑ

!!!!!!!!!! !!!!!!! !!!

“The more I learn, the more I realize I don’t know. The more I

realize I don’t know, the more I want to learn.”―Albert Einstein Abstract 本稿では，長記憶性をもつ時系列に対し，推定に関して最近 の研究を概観する。特に，長記憶性をもつ時系列で典型的な ARFIMAモデルについて，パラメトリックな推定法について 述べる。これらの推定法に対してシミュレーションでパフォー マンスを比較し，日本のファイナンス時系列データに応用する。 １ Introduction 長記憶性をもつ時系列に関しては多くの研究が行われている。長記憶 性を持つ時系列について概観したものには，Beran（１９９４），Robinson （２００３），Doukhan et. al.（２００３），Teyssiere and Kirman（２００７），Palma （２００７）がある。邦文文献としては，矢島（２００３）がある。

長記憶性については金融時系列を含む多くの分野で応用されている。 特に，金融時系列における実証的な結果としては，Ding et. al.（１９９３） によって指摘された通り，金融時系列のリターンそのままでは自己相関 を計算すると相関がないように見えるのに対し，リターンの絶対値およ び２乗の値の自己相関では，ラグを長くとるにつれての減少度合いが非

論 説

長記憶性をもつ時系列について

―ARFIMAモデルに対するパラメトリック推定

西 埜 晴 久

（３４７） １３９

常にゆっくりとしている。このことは，リターンの絶対値および２乗の 値が長記憶性を持つことを示しているといえる。 他方，長記憶性はHosking（１９８１）が示しているように，Fractional Difference（実数差分）と関連がある。Fractional Differenceを使って 表現される長記憶性をもつ典型的な時系列モデルとしては，ARFIMA （p，d，q）モデルがある。ARFIMA（p，d，q）モデルはBox＝Jenkins 法で用いられるARIMA（p，d，q）の階差パラメータd を整数値以外に も拡張したものである。本稿では，ARFIMA（p，d，q）モデルのパラ メトリックな推定法について解説することを目的とする。さらには，シ ミュレーションによって各推定法のパフォーマンスの比較を行う。 そして，取り上げた推定法を日本のファイナンス・データに実際に応 用する。本稿では，日経株価平均のrealized volatility（実現ボラティリ ティ）に対してARFIMAモデルを推定することを考える。Ding et. al. （１９９３）によれば，リターンそのものは無相関であるけれども，リター ンの２乗値および絶対値は非常にゆるやかな自己相関があることを示し ていた。このことは，リターンの分散を表している実現ボラティリティ にも長記憶性があることを示している。特に，渡部（２００７）では実現ボ ラティリティおよびその対数値に長記憶性があることを示していた。本 稿では，ARFIMAモデルの推定法を実現ボラティリティの対数値に対 して応用することを考える。 本稿は以下の内容で構成される。次の節では長記憶性とARFIMAモ デルについて説明する。３節では推定法について展望する。４節では推 定法に対しシミュレーションを行い，次に日経平均株価指数の実現ボラ ティリティのデータにそれらの推定法を応用する。 ２ 長記憶性とARFIMAモデル 定常時系列｛yt，t＝１，２，．．．，｝は自己共分散関数をフーリエ変換 １４０ （３４８）

することにより以下のスペクトル密度関数を持つ。 f（λ）＝ １_２π Σ∞ s＝−∞e isλ_{γ（s）． ¸} ただし，γ（s）＝Cov（yt，yt−s）は自己共分散関数である。 長記憶性は以下のようにスペクトル密度関数によって定義される。 f（λ）∼Cλ−２d_{，λ→＋０ （０＜d ＜１} ２）． ¹ このことはスペクトル密度は原点の近くで発散していることを示してい る。このことは同時に ２πf（０）＝ Σ∞ s＝−∞γ（s）＝∞， º を示しており，自己共分散の和が発散していることを示している。長記 憶性の定義として自己共分散の和が発散しているものを採ることもある。 スペクトル密度関数にせよ，自己共分散にせよ，定常過程においてのみ 定義できるものであるので，長記憶性は厳密には定常過程において定義 される概念であることに注意されたい。 次に，Fractional DifferenceをHosking（１９８１）に従って説明しよう。 まず， （１−L）d_y t＝εt，−１／２＜d＜１／２， » と書く。ただし，Lはラグ・オペレータであり，Lyt＝yt−１となる。また， ｛εt｝∼WN （０，σ２）は平均０分散σ２のホワイトノイズ過程である。し たがってMA表現を行えば， yt＝（１−L）−dεt， ¼ である。さらに二項展開より （１−L）−d_＝Σ∞ j＝０ Γ（j＋d ） Γ（j＋１）Γ（d ）Lj． ½ （３４９） １４１

である。ただし，Γ（x）＝

∫

∞_０zx−１_e−z_{dzはガンマ関数である。} したがって，｛yt｝のMA表現は以下のようになる： yt＝ ∞ Σ j＝０ Γ（j＋d ） Γ（j＋１）Γ（d ）εt−j． ¾ Fractional Differenceを用いることによって，以下のARFIMA（p，d， q）（autoregressive fractionally integrated moving average）モデルを 考えることができる。 φ（L）（１−L）d_{yt＝θ（L）εt，t＝１，２，．．．． ¿} ただし，｛εt｝はホワイトノイズプロセスWN （０，σ２）であり，φ（z）＝ ０とθ（z）＝０の全ての根は単位円外にあるものとする。 なお，０＜d ＜１_{２の時，｛}yt｝は定常な長記憶過程であり，他方d ＞_―１ ２ の時，｛yt｝は非定常な実数和分過程である。 ARFIMA（p，d，q）モデルは以下のスペクトル密度を持つ。 f（λ）＝ σ_{２π｜１−}２ eiλ_｜−２d｜θ（eiλ）｜２

｜φ（eiλ_）｜２． À

ここで

｜１−eiλ_｜−２d_{＝｜２sin（λ／２）｜}−２d_∼λ−２d_{，λ→＋０，}

（０＜d ＜１_２）， Á であることに注意すれば，ARFIMA（p，d，q）モデルが長記憶性を持 つことが分かる。なぜなら， f（λ）∼ σ_２π２ ｜θ_｜φ（１）_（１）｜_｜２２λ−２d，λ→＋０． Â となり，原点でのスペクトルが発散するからである。 例として，単純なARFIMA（１，d，１）モデルを考えてみよう。 １４２ （３５０）

（１−φL）（１−L）d_y

t＝（１−θL）εt， Ã

このときのARFIMA（１，d，１）モデルのスペクトル密度関数が原点 で発散することは以下のようにして確かめられる。

f（λ）＝ σ_２π｜２２ sin（λ_２）｜−２d｜１−θeiλ｜２

｜１−φeiλ_｜２ Ä ∼ σ_２π２ ｜１−θ｜_{｜１−φ｜}２２λ−２d，asλ→＋０． Å なお，スペクトル密度を推定するための統計量としてピリオドグラム が良く知られている。ピリオドグラムは以下のように定義される。 I（λ）＝ １T ２πT

｜

T Σ j＝１yje ijλ

｜

２_{． Æ} スペクトル密度を推定するピリオドグラムを使ったカーネル推定法とし て f ＾（λj）＝ h Σ m＝−h

［

h＋１−｜m｜２ （h＋１）２

］

I（λT j＋m）， Ç がある。ただし，周波数はλj＝２πj／T，（j＝１，２，．．．，T ）ととるも のとする。 ３ 推定法 本節ではARFIMA（p，d，q）モデルの推定法について説明する。は じめに，正規分布を仮定した時間領域の最尤法をSowell（１９９２）にした がい説明する。しかしながら，本推定法は計算にあたり難点がある。つ まり，標本サイズが大きくなるにつれて，自己共分散行列の逆行列を計 算することに多くの時間がかかってしまうのである。 こうした難点を回避するために時間領域の尤度を近似する方法がいく つか提案されている。良く使われているものとして偏自己相関を用いて， 一定の長さをとったあとは近似式を使うHaslett and Raftery（１９８９）が

ある。この方法は統計ソフトＲのパッケージにも使われているものであ る。その他に有力な方法として，AR近似を使ったBeran（１９９５）によ る推定法がある。

さらには，周波数領域の観点からはWhittle推定を使うことができる。 Whittle推定量に関しては，すでにFox and Taqqu（１９８６），Dahlhaus （１９８９）により長記憶性をもつ時系列においても，漸近正規性など漸近 的に良い性質を持つことが確かめられている。つまり，時間領域の厳密 な最尤推定量とWhittle推定量が漸近的には同等であることが分かって いる。

他には，d だけを推定するセミパラメトリックな方法が広く用いられ ている。例えば，Geweke et. al.（１９８３）によるGPH推定法およびRobin-son（１９９５）によるlocal Whittle推定法が代表的なものである。しかし， 本稿ではセミパラメトリックな推定法については扱わない。

３．１ 時間領域における最尤推定

ここでは，Sowell（１９９２）による厳密な最尤推定について述べる。ま

ず，y＝（y１，y２，．．．，yT）′が正規分布であると仮定する。すると，（y１，

y２，．．．，yT）の同時密度は，

f（y１，y２，．．．，yT｜d，φ，θ，σ２）＝

（２π）−T／２_｜σ２_Σ｜−１／２_exp

｛

_{− １}

２σ２y′Σ−１y

｝

， È

となる。よって，対数尤度は

logL（d，φ，θ，σ２_）＝−T

２log２π−１２log｜σ２Σ｜− １２σ２y′Σ−１y， É

となる。ただし，自己共分散行列は

σ２_{Σi,j＝γ（｜i−j｜） Ê}

σ２_Σ＝ ! % % % % % % % % % % % # γ（０） γ（１） … γ（T −１） γ（１） γ（０） … … … … … … γ（T −１） …… γ（０） " % % % % % % % % % % % $ となっている。 Sowell（１９９２）は¿式のARFIMA（p，d，q）モデルの自己共分散γ（s） を以下のように示している。 γ（s）＝σ２_Σq l ＝−q p Σ j＝１ψ（l ）ζjC（d，p＋l −s，ρj）． Ë なお，ψ（l ），ρj，ζj，C（d，h，ρ）は以下のように定義されている。 ψ（l ）＝min［q，q−l ］Σ s＝max［０，l ］θsθs−l， Ì φ（x）＝Πp j＝１（１−ρjx）， Í ζj＝

［

ρj p Π i＝１（１−ρiρj）Πm≠j（ρj−ρm）

］

−１ ， Î C（d，h，ρ）＝_{Γ（１−d ＋h）}Γ（１−２d ）_{Γ（１−d ）Γ}Γ（d ＋h）_{（d ）} ×［ρ２p ２F（d ＋h，１；１−d ＋h；ρ）１ −２F１（d −h，１；１−d −h；ρ）−１］． Ï ただし，２F（a，b；c；x）は超幾何関数であり，以下の式で定義される。１ ２F（a，b；c；x）＝１ ∞ Σ n＝０ （a）（b）n n （c）n xn n！ ＝１＋_１！cab x＋a（a＋１）b（b＋１）_{２！c（c＋１）} x２ ＋a（a＋１）_{３！c（c＋１）}（a＋２）b（b＋１）_{（c＋２）}（b＋２）x３_{＋…． Ð} （３５３） １４５

なお，（ ）nはPochhammerの記号と呼ばれ，（a）n＝a（a＋１）…（a＋n−

１）と定義される。また，超幾何関数の性質についてはAbramowitz and Stegun（１９６５）あるいはGradshteyn and Ryzhik（２０００）を参照の こと。

本推定法の難点は，Σの逆行列を計算するためにはCholesky分解を用 いることになるので計算量が非常に多くなり，標本サイズが大きくなる と計算に時間がかかることである。

３．２ Haslett-Raftery法

Haslett-Raftery法は，Haslett and Raftery（１９８９）によって提案され たものである。本推定法では，偏自己相関を使い全ての過去の値から一 期先予測値を計算する。ただし，計算量が増えることを防ぐため，用い る過去の値は一定の長さに制限する。実用上はM＝１００程度で良いよう である。 まず，以下の簡単なARFIMA（０，d，０）モデルを例としてHaslett-Raftery法を説明しよう。次のARFIMA（０，d，０）モデル： （１−L）d_y t＝εt，t＝１，２，．．．，T， Ñ を考える。ただし，εtはN（０，σ２ε）の正規変数とする。このARFIMA （０，d，０）にしたがうytに対し，一期先予測を考えると y ＾t＝ t−１ Σ j＝１φt,jyt−j， Ò となる。この時の予測誤差分散は， vt＝Var（yt−y＾t）＝σ２y t−１ Π j＝１（１−φ ２_{j,j）， Ó} となる。ただし， σ２_{y＝Var（yt）＝Γ}（１−２d ） Γ（１−d ）２σ２ε． Ô １４６ （３５４）

である。そこで，標準化された予測誤差分散rtは， rt＝ vt σ２ ε＝Γ （１−２d ） Γ（１−d ）２ t−１ Π j＝１（１−φ ２_{j,j）， Õ} r１＝Γ（１−２d ）_{Γ（１−d ）}２， Ö rt＝（１−φ２t−１，t−１）rt−１． × となる。 近似的な対数尤度は以下の通りである。 L＝−_２Tlog２π−_２Tlogσ＾２ ε−１ ２ T Σ t＝１logrt− T ２． Ø ただし， σ ＾２ ε＝１_T T Σ t＝１ （yt−y＾t）２ rt ． Ù となる。 あるいは以下のl を最大化すれば良い。 l ＝−１_２log

（

１_T ΣT t＝１ （yt−y＾t）２ rt

）

． Ú 次に，Òで用いる最良線形予測の係数｛φt,j｝を説明する。｛φt,j｝は偏 自己相関でもあり， φt,j＝−tCjΓ（j−d ）Γ（t−d −j＋１） Γ（−d ）Γ（t−d ＋１） ＝−Πj k＝１

（

（t−k＋１）（k−d −１） k（t−d −k＋１）

）

，j＝１，２，．．．，t． Û と求められる。｛φt,j；t＝１，．．．，T ，j＝１，．．．，M ｝を具体的に求 めるには，以下の再帰式を使って計算することになる。 φt，１＝_t−d，td Ü （３５５） １４７

φt,j＝φt,j−１

｛

（t−j＋１）（j−１−d ） j（t−d −j＋１）

｝

，j＝２，．．．t． Ý また， φj,j＝ d j−d， φt,t＝ d t−d， Þ である。 Ò_{の最後の項は以下のように近似する。} t−１ Σ j＝１φtjyt−j∼∼ M Σ j＝１φtjyt−j− t−１ Σ j＝M ＋１πjyt−j ß ∼ ∼ΣM j＝１φtjyt−j− MπM d

［

１−

（

M t

）

d

］

y―M＋１，t−１−M， à ただし。y―M ＋１，t−１−M＝_{t−１−２MΣ}１ t−１−Mj＝M＋１yjとおいており，打ち切り数 はM＝１００が望ましいとされている。 ß_{の近似式で用いられる係数π}jは，二項展開 （１−L）d_＝Σ∞ k＝０πjL j _á で得られるものである。そして，πjは以下の近似式を用いて計算する。 πj＝ Γ（j−d ） Γ（j＋１）Γ（−d ）∼∼ j−d−１ Γ（−d ） as j→∞， â πM＝ Γ（M −d ） Γ（M ＋１）Γ（−d ）∼∼ M−d−１ Γ（−d ） M is large． ã さらに，àの近似式は以下のように求められる。 t−１ Σ j＝M ＋１πj∼∼ t−１ Σ j＝M ＋１ j−d−１ Γ（−d ）∼∼Γ（−d ）１

∫

t M u −d−１_du ＝_{Γ（−d ）}１_{（−d ）}［t−d_−M−d_］ １４８ （３５６）

＝_{Γ（−d ）}M−d_{（d ）}

［

１−

（

M_t

）

d

］

∼ ∼Mπ_dM

［

１−

（

M_t

）

d

］

． ä 次に，ARFIMA（１，d，１）の場合のHaslett-Raftery法を説明しよ う。 ARFIMA（１，d，１）： （１−φL）（１−L）d_y t＝（１−θL）εt，t＝１，２，．．．，T． å ただし，εtはN （０，σ２εの正規変数である。ytの一期先予測は以下のよ うに計算され，そして近似される。 y ＾t＝（１−θL）−１（１−φL） t−１ Σ k＝１φtkyt−k ∼ ∼（１−θL）−１_{（１−φL）}

｛

_ΣM k＝１φtkyt−k− MπM d

［

１−

（

M t

）

d

］

y―M＋１，t−１−M

｝

＝（１−θL）−１min（M，t−２）_Σ k＝１ φ（ytk t−k−φyt−１−k） −

（

１−φ_１−θ

）

Mπ_dM

［

１−

（

M_t

）

d

］

y―M＋１，t−１−M ＝min（M，t−１）Σ j＝１ φtj min（M，t−j−２） Σ k＝０ θ k ｛yt−j−k−φyt−j−k−１｝ −

（

１−φ_１−θ

）

Mπ_dM

［

１−

（

M_t

）

d

］

y―M＋１，t−１−M． æ そして予測誤差分散は vt＝Var（yt−y＾t）＝σ２yκ t−１ Π j＝１（１−φjj） ２_{， ç} となる。ただし， σ２_{y＝Var（yt）＝Γ}（１−２d ） Γ（１−d ）２σ２ε， è （３５７） １４９

κ＝_１＋θ１−φ２_{−２φθ．}２ é である。なおκはARMA（１，１）モデルのイノベーションの分散と ARMA（１，１）自体の分散との比である。 よって標準化された予測誤差分散rtは以下のようになる。 rt＝vt σ２ ε＝κΓ （１−２d ） Γ（１−d ）２ t−１ Π j＝１（１−φjj） ２_{， ê} r１＝

（

１−φ ２ １＋θ２_−２φθ

）

Γ（１−２d ）_{Γ（１−d ）}２ ， ë rt＝（１−φ２t−１，t−１）rt−１． ì したがって， σ ＾２ ε＝１_T T Σ t＝１ （yt−y＾t）２ rt ． í を最小化することで，d，φ，θを推定できる。あるいは， l ＝−logσ＾２ ε． î を最大化して推定する。 最後に一般のARFIMA（p，d，q）の場合のHaslett-Raftery法を説明 しよう。 ARFIMA（p，d，q）モデル： φ（L）（１−L）d_y t＝θ（L）εt，t＝１，２，．．．，T． ï を考えると，ytの一期先予測は次のように計算される。 y ＾t＝θ（L）−１φ（L） t−１ Σ k＝１φtkyt−k． ð よって，これをM で打ち切って近似してy＾tを求めてやればよい。そして 予測誤差分散は vt＝Var（yt−y＾t）＝σ２yκ t−１ Π j＝１（１−φjj） ２_{， ñ} １５０ （３５８）

となる。ただし， σ２_{y＝Var（yt）＝Γ}（１−２d ） Γ（１−d ）２σ２ε， ò であり，κはARMA（p，q）モデルのイノベーションの分散とARMA （p，q）自体の分散との比である。よって標準化された予測誤差分散rt は以下のようになる。 rt＝ vt σ２ ε＝κΓ （１−２d ） Γ（１−d ）２ t−１ Π j＝１（１−φjj） ó したがって，ytの一期先予測y＾tおよび予測誤差分散rtを使って， σ ＾２ ε＝１_T T Σ t＝１ （yt−y＾t）２ rt ． ô が計算される。このσ＾２ εを最小化することでパラメータを推定できる。 あるいは， l ＝−logσ＾２ ε． õ を最大化して推定する。 ３．３ BeranのAR近似法 Beran（１９９５）はARFIMAモデルをAR近似することで推定する方法 を提案した。まずは，ARFIMA（０，d，０）を例として推定法を説明 する。 次のARFIMA（０，d，０）モデル： （１−L）d_y t＝εt，t＝１，２，．．．，T． ö を考える。ただし，εtはN （０，σ２ε）の正規変数である。 ここで二項展開より （１−L）d_＝Σ∞ k＝０πkL k _÷ であり，その係数は （３５９） １５１

πk＝ （k−d −）！ k！（−d −１）！＝Γ（k＋１）ΓΓ（k−d ）（−d ） ＝Πk j＝１ j−１−d j ． ø と求められる。よって， π０＝１ ù π１＝−d ú π２＝（１−d ）_２（−d ） û πk−１＝ k−１ Π j＝１ j−１−d j ü πk＝

（

k−１−d_k

）

πk−１． @ とすることで再帰的に計算することができる。 したがって，öは εt＝ ∞ Σ k＝０πkyt−k，t＝１，２，．．．，T， A （ただしπ０＝１）と計算することができる。 この時の残差は et＝ t−１ Σ k＝０πkyt−k＝yt＋ t−１ Σ k＝１πkyt−k，t＝１，２，．．．，T． B と計算することができる。すなわち， e１＝y１ C e２＝y２＋π１y１ D e３＝y３＋π１y２＋＋π２y１ E et＝yt＋ t−１ Σ k＝１πkyt−k F １５２ （３６０）

eT＝yT＋ T −１ Σ k＝１πkyT−k． G と計算される。もし，期待値E（yt）が未知ならば， et＝yt−y―＋ t−１ Σ k＝１π（yk t−k−y―），t＝１，２，．．．，T． H と平均を引いて計算すればよい。 残差｛et｝によって近似対数尤度が計算できる。したがって Lberan＝−T ２log２π− T ２logσ２ε−１_２ T Σ t＝２

｛

et σε

｝

２ I を最大化すれば良い。 もし，以下のARFIMA（１，d，０）： （１−φL）（１−L）d_y t＝εt，t＝１，２，．．．，T． J を推定するならば， （１−φL）（１−L）d_{＝（１−φL）Σ}∞ k＝０πkL k ＝Σ∞ k＝０πkL k_−Σ∞ k＝０φπkL k＋１ ＝１＋Σ∞ k＝０（πk−φπk−１）L k _K の式を用いれば良い。つまり，π∼ k＝πk−φπk−１とおけば， π∼ k＝（k−d −１）！ k！（−d −１）！−φ（k−１）！（−d −１）！（k−d −２）！ ＝

（

k−d −１_k −φ

）

_{（k−１）！（−d −１）！}（k−d −２）！ ＝

（

１−d ＋１_{k −φ}

）

Γ_Γ（k−d −１）_{（k）Γ（−d ）} （３６１） １５３

＝

（

１−d ＋１_{k −φ}

）

k−１Π j＝１ j−１−d j ， L となるので，残差を et＝yt＋ t−１ Σ k＝１π ∼ kyt−k，t＝１，２，．．．，T． M と計算する。後は，Iによって計算すれば良い。 次に，ARFIMA（０，d，１）： （１−L）d_y t＝（１−θL）εt，t＝１，２，．．．，T， N の場合を説明する。 （１−θL）−１_{（１−L）}d_＝

（

_Σ∞ j＝０θ j_Lj

）（

_Σ∞ k＝０πkL k

）

＝Σ∞ k＝０

｛

k Σ j＝０πjθ k−j

｝

_Lk ＝１＋Σ∞ k＝１

｛

k Σ j＝０πjθ k−j

｝

_Lk _O であるから，π∼k＝Σkj＝０πk−jとおくことで，残差｛et｝が et＝yt＋ t−１ Σ k＝１π ∼ kyt−k，t＝１，２，．．．，T． P と計算できる。後は，Iを使えば良い。 そして，ARFIMA（１，d，１）： （１−φL）（１−L）d_y t＝（１−θL）εt，t＝１，２，．．．，T Q の場合には， （１−θL）−１_{（１−φL）（１−L）}d ＝

（

Σ∞ j＝０θ j_Lj

）

（１−φL）

（

Σ∞ k＝０π０L k

）

＝

（

Σ∞ j＝０θ j_Lj

）［

∞ Σ k＝０πkL k_−φΣ∞ k＝０πkL k＋１

］

１５４ （３６２）

＝Σ∞ k＝０

｛

k Σ j＝０πjθ k−j

｝

_Lk_−φΣ∞ k＝０

｛

k Σ j＝０πjθ k−j

｝

_Lk＋１ ＝１＋Σ∞ k＝１

｛

k Σ j＝０πjθ k−j_−φk−１_Σ j＝０πjθ k−１−j

｝

_Lk_{． R} を使って計算する。つまり， π∼k＝ k Σ j＝０πjθ k−j_−φk−１_Σ j＝０πjθ k−１−j _S とおくことで，残差｛et｝ et＝yt＋ t−１ Σ k＝１π ∼ kyt−k，t＝１，２，．．．，T． T を計算して，Iを使えば良い。 最後に一般のARFIMA（p，d，q）の場合を説明しよう。ARFIMA （p，d，q）モデル： φ（L）（１−L）d_y t＝θ（L）εt，t＝１，２，．．．，T． U を考える。この場合も， θ（L）−１_{φ（L）（１−L）}d_＝Σ∞ j＝０π ∼ jLj V の展開式を使って，残差｛et｝ et＝yt＋ t−１ Σ k＝１π ∼ kyt−k，t＝１，２，．．．，T． W を求める。そして， Lberan＝−T ２log２π− T ２logσ２ε−１_２ T Σ t＝２

｛

et σε

｝

２ X を最大化すれば良い。 ３．４ Whittle推定（周波数領域） Whittle推定は時間領域の尤度を周波数領域で表現することで求めら （３６３） １５５

れる。詳しくはBrockwell and Davis（１９９１）あるいは谷口（２００３）を 参考にしていただきたい。 まず，ARFIMA（p，d，q）のパラメータをτ＝（φ，θ，d，σ２_{）′とし，} 周波数をλj＝２π_{T ，}j j＝１，．．．，T −１ととる。そして，I（λT j）をピリオ ドグラムとすると，Whittle推定量τ∼W_は， LW_{（τ）＝ １} ２T T −１ Σ j＝１

｛

logf（λτ j）＋ I（λT j） fτ（λj）

｝

， Y を最小化することで得られる。 Dahlhaus（１９８９）は長記憶性を持つ時系列の場合にも一致性と漸近 正規性を示している。そして，Whittle推定は以下の式より時間領域の 最尤推定の近似になっていることが分かる。 −TLW_{（τ）＝−１} ２ T −１ Σ j＝１logfτ（λj）−１２ T −１ Σ j＝１ I（λT j） fτ（λj） Z ∼ ∼const−１_２log｜σ２_{Σ｜− １} ２σ２y′Σ−１y． [

Chueng and Diebold（１９９４）は３．２節で述べた時間領域の最尤法と Whittle推定のパフォーマンスをシミュレーションによって比較してい る。そこでは，標本サイズがある程度大きい場合にWhittle推定は時間 領域の最尤法と同程度のパフォーマンスであるとしている。また平均が 未知の場合には，平均を推定する必要のないWhittle推定が有用であろ うとしている。 最後にARFIMA（p，d，q）をWhittle法で推定する際のreparametri-zationについて述べておこう。ARFIMA（p，d，q）のスペクトル密度 は

f（λ）＝ σ_{２π｜１−}２ eiλ_｜−２d｜θ（eiλ）｜２

｜φ（eiλ_）｜２． \

であるから，θ（z）とφ（z）を θ（z）＝Πq j＝１（１−θjz）and φ（z）＝ p Π i＝１（１−φiz）， ] と分解できる形に表現してやると，｜φi｜＜１および｜θj｜＜１が定常性 および反転可能性を示すので扱いやすくなる。そこで ｜θ（eiλ_）｜２_＝Πq j＝１｜１−θje

iλ_）｜２_and｜φ（eiλ_）｜２_＝Πp

i＝１｜１−φie iλ_｜２ _^ を使ってスペクトル密度を計算してやればよいことになる。 ４ シミュレーションとファイナンス・データへの応用 本節では前節でとりあげた推定法について，まずシミュレーションで そのパフォーマンスを確認する。次にファイナンス・データである日経 平均株価指数のrealized volatilityに対して，実際に応用する。 ４．１ シミュレーション 本 節 で は 前 節 で 取 り 上 げ た４つ の 推 定 法：exact MLE，Haslett-Raftery法，BeranのAR近似法，Whittle推定について，モンテカルロ実 験によってそのパフォーマンスを比較する。本実験では，ARFIMA（０， d，０）モデルにしたがう時系列データを発生させる。d の値はd ＝０．２ およびd ＝０．４とする。d ＝０．４の方が長記憶性の程度が強いといえよう。 また，標本サイズはT ＝２００，T ＝５００，T ＝１，０００の三つを用いた。複数 の標本サイズを比較することで，標本が大きくなるときの漸近的な特性 を見ることができるだろう。また，実際の値としては，推定値のバイア スおよび平均２乗誤差の平方根（Root Mean square error，RMSE）を 計算した。バイアスおよびRMSEが小さいものが良い推定量となる。

Table１にモンテカルロ実験の結果がまとめてある。左の列からexact ML，Haslett-Raftery法，Beran（１９９５）のAR近似法，Whittle推定法の

４つである。なお，Haslett-Raftery法ではM ＝１００で打ち切って計算し， Beran（１９９５）のAR近似法ではm＝５０で打ち切って計算している。上段 にはバイアスが下段にはRMSEが記されている． まず，exact MLはT ＝２００では最適化に失敗してしまうことがあった。 さらに，T ＝１，０００の時のexact MLは非常に計算時間がかかってしまう。 exact MLは自己共分散行列の逆行列を得ることに時間がかかるのが難 点であるためであろう。しかしながら，RMSEで判断すると，exact MLが一番良い。その次はHaslett-Raftery法が全体的には良いであろう。 しかし，d ＝０．４と非定常な領域に近づく時およびT ＝１，０００のように標 本サイズが大きいときには，Whittle推定法のRMSEが小さく，良いパ

Table １ Estimates ofd of ARFIMA（０,d,０）

Exact H-R AR approx Whittle T ＝２００ d ＝０．２ −０．００７４０６１ −０．００２４６９６ −０．００３６６９９ −０．００４３４９８ ０．０５５１３２ ０．０５７７２３ ０．０５８２４３ ０．０６１７８７ d ＝０．４ −０．０１１４５７ ０．０１３７７５ ０．０１５１２３ −０．０００８９２６９ ０．０４９１９９ ０．０５９６３７ ０．０６０４８８ ０．０６００５０ T ＝５００ d ＝０．２ −０．００２７０６９ −０．０００２４１０２ ０．０００６３０１６ −０．００１７８８９ ０．０３５２９５ ０．０３５８７８ ０．０３６５６１ ０．０３７７３７ d ＝０．４ −０．００３５１５８ ０．００９０１７５ ０．０１９３０３ ０．００１１６６３ ０．０３２７９７ ０．０３８５１５ ０．０４６０８５ ０．０３７９６６ T ＝１，０００ d ＝０．２ −０．００１８７３３ −０．０００４１７４２ ０．００１４４９３ −０．００１６４３２ ０．０２５４０５ ０．０２５６８３ ０．０２６５１６ ０．０２６４６３ d ＝０．４ −０．００１８７３０ ０．００５６４８７ ０．０２０５３２ ０．０００５９１２８ ０．０２４２３７ ０．０２６９８６ ０．０３９８３６ ０．０２６５９６ upper: Bias, lower: RMSE

フォーマンスを示している。また，Whittle推定法のバイアスも全体的 に小さいといえるだろう。BeranのAR近似法はd ＝０．２０の時および標本 サイズが小さいときには，Haslett-Raftery法についでRMSEが小さい。 一方，BeranのAR近似法標本サイズが大きくなってもHaslett-Raftery 法あるいはWhittle推定法ほどにはRMSEが小さくはならない。 ４．２ ファイナンス・データへの応用 これまで述べてきた推定法を日本のファイナンス時系列にデータに応 用する。ここで利用する日経平均株価指数は日本の株式市場で最も典型 的な株価指数である。本稿では日経平均の５秒毎のデータからrealized volatility（実現ボラティリティ）を計算して，その対数値に対してAR-FIMAモデルを推定することを目指す。 まず，日次realized volatility，RVtは以下のように定義されるもので ある。 RVt＝c n Σ i＝１r ２_{i,t， １０１½} c＝ΣTt＝１（Rt−R―）２ ΣT t＝１Σni＝１r２i,t ½ １０２ なお，Rtは日次リターンであり，R―は日次リターンの平均である。そし て，ri,tはt日における５秒毎の対数収益率であり，cは½１０２で定義される調

整定数である。（例えば，Hansen and Lunde（２００５），渡部（２００７）を 参照のこと）。

なお，本稿で用いる日 経 平 均 株 価 指 数 の デ ー タ は，標 本 期 間 が

２００１．７．２―２００８．６．３０であり，標本サイズはT ＝１，７０７，調整定数は

c＝２．０２２である。

Table２はrealized volatility RVt，realized volatilityの対数値 log（RVt），

日次リターンRt，日次リターンの２乗値R２tの要約統計量を示している。

RtのLjung-Box統計量は小さく，そのp値は０．５９３４である。このことは，

Rt自体には系列相関がないという帰無仮説を棄却することはできない。

一方，log（RVt）のLjung-Box統計量は非常に大きく，高い系列相関が

あることが考えられる．

Figure１は日経平均株価指数から計算されたrealized volatilityをグラ Table２ Summary statistics, Nikkei

RVt log（RVt） Rt Rt２ mean １．９７１８０ ０．３９２９７ ０．００２２７ １．９７１８１ stdev １．６１８１７ ０．７８０７８ １．４０４６２ ３．６９６３３ skewness ２．６８２８９ −０．２３１４９ −０．３１１６６ ５．２１４２１ kurtosis １６．６６４３６ ２．８５１８４ ４．５０８７７ ４６．６８１７１ LB（１０） ４，５６５．５４ ７，４５７．２８ ８．３６３５ ２５６．４５５ Rt: Daily Return.

LB（１０）:Ljung-Box test statistic with lag１０.

Figure１ RV of Nikkei２２５,２July２００１―３０June２００８

フにしたものである。一方，Figure２はrealized volatilitiyとその対数値 の自己相関をグラフにしたものである。このグラフからは，realized volatilityの対数値の自己相関の減衰が非常に遅いことが分かる。つまり， realized volatilityの対数値系列は長記憶性を持っていると考えられるの である。

次にFigure３でrealized volatilityの対数値系列 log（RVt）とその推定

されたスペクトル密度をグラフにしている。この推定されたスペクトル 密度は原点で発散している。このことは log（RVt）が長記憶時系列か非 定常時系列であることを示している。すなわち，log（RVt）に対して何 らかの階差をとることが必要と考えられるのである。 一方，Figure４はrealized volatilityの対数の１階の階差系列とその推 定されたスペクトル密度をグラフにしている。この推定されたスペクト ル密度は，１階の階差をとることは過剰階差であることを示している。 したがって，Figure３，４を合わせて考えれば，realized volatilityの対

Figure２ Autocorrelations of RV and log（RV）,Nikkei

Figure３ RVの対数値系列とその推定されたスペクトル密度

Figure４ RVの対数値の階差系列とその推定されたスペクトル密度

数値に対しては実数差分をとることが望ましいと考えられるのである。 Table３にはHaslett-Raftery法，AR近 似 法，Whittleの３つ の 推 定 法 で推定した結果を表にしてある。RFIMA（０，d，０）以外のモデルで は，数値的最適化がうまくいかなかったので，ARFIMA（０，d，０） の推定結果のみを表に記してある。またexact MLでは標本サイズが大 きいためARFIMA（０，d，０）の推定の計算ができなかった。Table３ には推定値と標準誤差（S.E.）そしてt値が記されている。d の推定値は d ＝０．４５８―０．４８４であり，極めて０．５に近い値である。また，t値が大き いので十分に有意と判断できるだろう。 なお，３つの推定法の中で中間の値をとったHaslett-Raftery法での推 定値d ＝０．４６０によって，realized volatilityの対数値系列に対してフィル ターをかける。フィルターをかけた後の推定されたスペクトル密度は Figure５に描かれている。このスペクトルは原点で発散していないので，

Table３ Estimates of ARFIMA（０,d,０）,Nikkei

Haslett-Raftery AR approximation Whittle estimates of d ０．４６０１３ ０．４８４１４ ０．４５８０５ S.E.s ０．０１１７７０ ０．０１６０３８ ０．０１２０９５ t-values ３９．０９２ ３０．１８７ ３７．８７１

Figure５ d＝０．４６０でフィルターをかけたあとの推定されたスペクトル密度

長記憶時系列ではなく，また，過剰階差でもないであろう。つまり，短 記憶定常時系列になっていると考えられる。 謝辞 本研究は，財団法人日本証券奨学財団および財団法人野村財団，そし て科学研究費補助金＃２１５３０１９４の助成を受けた。また，プログラムの一 部の作成において各務和彦氏の協力を得た。記して感謝申し上げる。も ちろん，本稿のありうべき誤りの責任は全て著者にある。 References ［１］ 谷口正信（２００３）「時系列解析入門」，竹内啓編『統計学の基礎―線形モデルからの出 発』岩波書店所収。 ［２］ 矢島美寛（２００３）「長期記憶をもつ時系列モデル」，竹内啓編『経済時系列の統計―そ の数理的方法』岩波書店所収。 ［３］ 渡部敏明（２００７）「Realized Volatility――サーベイと日本の株式市場への応用――」， 『経済研究』，５８，３５２―３７３。

［４］ Abramowitz, M. and Stegun, I.A.（１９６５）（eds）Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, Dover, New York.

［５］ Beran, J.（１９９４）Statistics for Long-Memory Processes, Chapman & Hall, London. ［６］ Beran, J.（１９９５）“Maximum Likelihood Estimation of the Differencing Parameter for

Invertible Short and Long Memory Autoregressive Integrated Moving Average Mod-els”,Journal of the Royal Statistical Society. Ser. B ５７６５９―６７２.

［７］ Brockwell, P.J. and Davis, R.A.（１９９１）Time Series: Theory and Methods, ２nd ed., Springer, New York.

［８］ Chan, N.H. and Palma, W.（１９９８）“State Space Modeling of Long-Memory Processes”, Annals of Statistics ２６７１９―７４０.

［９］ Cheung, Y.-W. and Diebold, F.X.（１９９４）“On maximum likelihood estimation of the differencing parameter of fractionally-integrated noise with unknown mean”,Journal of Econometrics ６２,３０１―３１６.

［１０］ Dahlhaus, R.（１９８９）“Efficient parameter estimation for self-similar processes”,

nals of Statistics １７,１７４９―１７６６.

［１１］ Ding, Z., Granger, C.W.J. and Engle, R.F.（１９９３）“A long memory property of stock market returns and a new model”,Journal of Empirical Finance １,８３―１０６.

［１２］ Doukhan, P., Oppenheim, G. and Taqqu, M.S.（eds.）（２００３）Long-Range Dependence: Theory and Applications, Birkhauser, Boston.

［１３］ Fox, R and Taqqu, M.S.（１９８６）“Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series”,Annals of Statistics １４,５１７―５３２. ［１４］ Geweke, J. and Poter-Hudak, S.（１９８３）“The estimation and application of long

mem-ory time series models”,Journal of Time Series Analysis ４,２２１―２３８.

［１５］ Gradshteyn, I.S. and Ryzhik, I.M.（１９９４）Tables of Integrals, Series, and Products, ５th ed., Academic Press, New York.

［１６］ Granger, C.W.J. and Joyeux, R.（１９８０）“An introduction to long memory time series model and fractional differencing”,Journal of Time Series Analysis １,１５―３９.

［１７］ Hansen, P.R. and Lunde, A.（２００５）“A Comparison of Volatility Models: Does Any-thing Beat a GARCH（１，１）？”,Journal of Applied Econometrics ２０,８７３―８８９.

［１８］ Haslett, J. and Raftery, A.E.（１９８９）“Space-Time Modelling with Long-Memory De-pendence: Assessing Ireland’s Wind Power Resource”,Journal of the Royal Statistical So-ciety. Ser. C ３８,１―５０.

［１９］ Hosking, R.M.J.（１９８１）“Fractional differencing”,Biometrika ６８,１６５―１７６. ［２０］ Palma, W.（２００７）Long-Memory Time Series: Theory and Methods, Wiley, New York. ［２１］ Robinson, P.M.（１９９５）“Gaussian semiparametric estimation of long range

depend-ence”,Annals of Statistics,２３,１６３０―１６６１.

［２２］ Robinson, P.M.（ed.）（２００３）Time Series with Long Memory, Oxford, Oxford.

［２３］ Sowell, F.（１９９２）“Maximum likelihood estimation of stationary univariate fraction-ally integrated time series models”,Journal of Econometrics,５３,１６５―１８８.

［２４］ Teyssiere, G. and Kirman, A.P.（eds.）（２００７）Long Memory in Economics, Springer, New York.

（２０１０年９月１３日受理）

Summary

On Time Series with Long Memory

―Parametric Estimation of ARFIMA Model

Haruhisa NISHINO

This paper surveys the recent research on estimation methods for long memory time series. It focuses on parametric estimations for the ARFIMA model, which is typical among time series models with long memory. The paper moreover conducts a simulation study to compare the performance of the estimation methods, which are also applied for Japanese financial time series data.