山形大学大学院社会文化システム研究科紀要

第15号（２０１8）31－44

モンテカルロフィルタの理論と応用――SCD モデルへの応用――
砂田（社会システム専攻

１

洋志

企業経営領域担当）
方法であるために計算負荷が高い。したがって，

はじめにコンピュータの性能が高くないと推定できない。

経済現象をモデル化する際に用いられる経済変このため今まで広く利用されてこなかったと考え数の中には観察可能な変数がある一方で，観察不られる。今後も経済理論やモデルが複雑化される可能な変数もある。観察不可能な変数は状態変数，

ものの，コンピュータの性能は向上し続ける。し状態変数を用いたモデルは状態空間モデルと呼ばたがって，複雑なモデルの推定を可能にするモンれる。景気のように，状態変数の中には観察可能

テカルロフィルタへの期待は今後，一層高まるとな経済変数に大きな影響を与えるものがある。そ考えられる。

のため，状態空間モデルが必要とされ，経済学を日本の先行研究を調べた限り，経済現象を状態始めとして多くの分野で利用されている。

空間モデルによって定式化した上でモンテカルロ

状態空間モデルは，観測可能な変数と状態変数フィルタを用いてパラメータや状態変数を推定しの関係を記述した観測方程式，状態変数の時系列た研究は多くないと考えられる。先行研究として，的な関係を記述した遷移方程式によって構成され金利モデルに適用した佐藤 ･ 高橋（2005），確率る。状態変数の推定は多くの分野で関心を持たれ

ボラティリティーモデルへ適用した矢野 ･ 佐藤ているが，当初は観測方程式と遷移方程式が線形（2007），NAIRU（インフレを加速させない失業である上，両方の方程式内に含まれる誤差項が正率）に適用した高部（2010），DSGE（動学的確

規分布に従う場合に限って，Kalman（1960）で率的一般均衡）モデルの推定に適用した矢野 ･ 飯

提案されたカルマンフィルタを用いて推定するこ

田 ･ 和合（2011）などを挙げることができる。

と が で き た。 統 計 的 な 手 法 が 進 化 し て，

本稿では，モンテカルロフィルタをデュレー

Kitagawa（1993，1996） と Gordon et.al.（1993）

ションの時系列的な分析へ応用する。デュレー

によって線形性と正規性が満たされない状態空間

ションとは，証券市場や商品市場の１日の取引のモデルであっても推定できる方法が提案された。

時間間隔であり，１日の中でパターンを持っていKitagawa（1993，1996）においてこの方法はモンる。寄り付きと大引けは取引が集中するので取引

テカルロフィルタ（monte carlo filter）と呼ばれ間隔が短くなるものの，昼間は取引が少なくなる

1

ている。

モンテカルロフィルタはさらに発展し，ので取引間隔が広がるのである。こうしたパター

Kitagawa（1998）において未知パラメータも状

ンは日中効果（intraday effect）と呼ばれている。

態変数と考えて推定する自己組織化状態空間モデ

時系列モデルを用いたデュレーションの研究としルが提案された。

て ， E n g l e = R u s s e l l （ 1997 ， 1998 ） の A C D

モンテカルロフィルタはカルマンフィルタより（Autoregressive Conditional Duration）モデルも一般的な推定方法であるが，数値計算的な推定を基礎として様々なモデルが提案されてきた。2
その中で本稿では発展的なモデルである

1

Kitagawa（1993，1996）ではモンテカルロフィルタという名称を用いているが、近年は粒子フィルタ（particle
filter）という名称も広く用いられている。

2

先行研究については Pacurar（2008）に詳しくまとめら
れている。

－31－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田洋志）

Bauwens=Veredas（2004）の SCD（Stochastic

そして，ベイズ統計学を用いて SCD モデルを推Conditional Duration） モ デ ル を 取 り 上 げ る。

定した砂田（2018）の結果と比較する。第6節でSCD モデルの先行研究としては，擬似最尤法では本稿の内容をまとめる。

推定する Bauwens=Veredas（2004）に加えて，

２

ベイズ統計学の立場から SCD モデルを推定する

モンテカルロフィルタ3

方法を示した Strickland 他（2006），レバレッジ

状態空間モデルは観測方程式と遷移方程式で構効果を導入した SCD モデルをスライスサンプ

成される。両方の方程式が線形の式で記述される

ラーによって推定する Men 他（2015，2016）など

と共に誤差項が正規分布に従う状態空間モデルでがある。日本の市場を対象としたデュレーション

あれば，線形正規状態空間モデルと呼ばれ，状態の研究に限定すると森保（2006）が挙げられる。

変数について解析的な解を得ることができる。し森保（2006）では Log － ACD モデルと ACM －かし，そうでなければ解析的に解を得ることがで

ACD モデルを用いて東京金先物市場が計量的にきないので，数値計算的に解を得ることとなる。

分析されている。杉浦 = 中妻（2016）では，株

線形性や正規性が仮定されない状態空間モデルは式市場のティックデータへ SCD モデルが適用さ一般状態空間モデルと呼ばれ，その推定にはモンれ，モンテカルロフィルタを用いてパラメータがテカルロフィルタが利用される。この節では，一推定されている。さらに，推定結果を用いて約定

般状態空間モデルの状態変数を数値計算的に推定時間の予測も行なわれている。本稿は杉浦 = 中する方法であるモンテカルロフィルタについて説

妻（2016）を大いに参考にしている。この他には明する。

幾つかの研究が行われているに過ぎない。
2.1

前述したとおり，モンテカルロフィルタは非線

一般状態空間モデル

形な状態空間モデルや誤差項が正規分布に従わな

モンテカルロフィルタで扱う一般状態空間モデ

い状態空間モデルのパラメータや状態変数の推定

ルは線形モデルとは限らないので，以下の式のよに利用できる有益な推定方法である。そこで，モ

うな関数形で記述する。関数形で記述すると，線

ンテカルロフィルタを紹介することが本稿の第1形モデルは特殊形として含まれる。

の目的である。そして，第2の目的は現実の経済

遷移方程式

モ デ ル へ の 応 用 で あ り， フ ァ イ ナ ン ス 分 野 の観測方程式

SCD モデルを取り上げて予測・フィルタリング・

た だ し，

（1）（2）

は 観 測 さ れ る 時 系 列 デ ー タ，

平滑化という3種類の過程を紹介することである。

は観測されない状態変数である。

本稿の構成は以下のとおりである。まず，第2観測方程式における誤差項であり，その確率密度である。

は節ではモンテカルロフィルタについて説明する。

関数はは遷移方程式におけるその際，自己組織化状態空間モデルの推定方法も誤差項であり，その確率密度関数は紹介する。第3節では応用例の SCD モデルを紹

以下の小節において一期先予測とフィルタリング

介する他，自己組織化状態空間モデルを用いて同の方法を説明する。

である。

モデルを推定する方法を説明する。第4節では
2.2

データの説明とデータを平滑化する方法について説明する。第5節では，東京商品取引所に上場さ

一期先予測

一期先予測を行なうには，れている金先物市場のティックデータを対象にンプリングして得た状態変数の粒子

SCD モデルを適用して，状態変数とパラメータ

3

からサ
に

この節の記述は北川（2005）の付録D.
1とD.2，樋口編（2011）の6.2節に大きく負っている。

をフィルタリングによって推定した結果を示す。
－32－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田加 え て，洋志）

からサンプリングして得た粒子
も

必

要

で

あ

る。

な



お，

は時点 t までの時系列デー

一方，は以下のように式変形できる。

タ で あ る。 こ れ ら の 粒 子 を（1）の 遷 移 方 程 式

に代入することによって，状態変数 xtの一期先予測の粒子（3）を得ることができる。

（4）

そこで，（3）を（4）へ代入することによって次式が得られる。

そ の 理 由 は，の従う確率分布

が以下のとおりに記述できるからである。

の場合にだけ（5）と な る か ら，
となる。
また，
ディラッ

クのデルタ関数の性質から
ただし，δ はディラックのデルタ関数である。

である。以上の点を（5）へ代入すると次式が得ら
れる。

2.3

フィルタリング

フィルタリングに必要な尤度の説明から始める。状態変数 xt の一期先予測分布
ら生成された粒子

かが xt と等しい場合，尤度はであり，これをと記述する。

測方程式

の逆関数

に変更した

を

は，観

の xt をへ代入することによっ

上式はて，以下の式で記述される。

がの各粒子に重みをつけたもので近似できることを示している。さらに，を満たす整数列を導入すると，上式は以下のように記述できる。
から生成される粒子

は，から直接生成されるのではなく，上述した尤度を利用して一期先予測分布から生成される粒子を以下のよ
（6）に従えば，の複製が（6）個ずつ含まれるような粒子を生成する，あるいは，にうにリサンプリングすることで得られる。具体的含まれる各粒子が尤度にはされるように m 個の粒子を復元抽出することによって，を生成できることになる。つまり，を用いてことによって
2.4

に従って選択すれば良い。その理由は以下のとおりである。

に比例する割合で抽出をリサンプリングするが得られるのである。

平滑化

この小節では3種類の平滑化の方法を紹介する。

を用いると予測分布は以下のように近平滑化には，一般的な意味での平滑化（固定区間似される。

平滑化）に加えて，固定点平滑化と固定ラグ平滑－33－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田化の2種類がある。まず，時点1から時点 n までの全データが与えられた上での時点おける状態変数の確率分布洋志）

リ サ ン プ リ ン グ に 用 い る 粒 子にはが固定区間平…と逐次的に粒子を生成してい滑化分布である。次に，時点1から時点 t までのくことで得られる。そして，

データが与えられた上での時点移方程式（1）へ代入することで得られる。

る状態変数の確率分布におけが固定点平滑化分

はグすることでタが与えられた上での時点 t における状態変数のきる理由は以下のとおりである。

が固定ラグ平滑化分布である。

時点 t －1までの観測値

樋口編（2011）によれば，状態変数ベクトルが高与えられた場合，状態変数次元な場合の平滑化分布に関する推測は，固定区

時確率密度

間平滑化分布以下のように更新される。

分布を遷をリサンプリン布である。最後に，時点1から時点 t+L までのデー
確率分布との代わりに固定ラグ平滑化を得ることがでに加えて観測値 がに関する同はを用いて実行されると記述されている。以下では，利用されることが少なかったと思われる固定点平滑化と固定ラグ平滑化の手順について説明する。
⑴

固定点平滑化
固定点平滑化によって得られる粒子
はの j 番目の実現値である。つまり，時点 t 期までの情報 を既，（3）と知とした場合の時点 t までの状態変数から，に 関 す る j 番 目 の 粒 子 で あ る。

である。したがって，時点 t までの情報 に基づは，時点 t －1までの情報

いて平滑化，つまり，確率分布を 与 え ら れ た 場 合 の 時 点 t －1ま で の 状 態 変 数
の同時分布

から得

られる粒子
る

粒

子

個 の t次 元 ベ ク ト ル の 粒 子をフィルタリングの場合と同じ重みから得らでリサンプリングすれば良いのである。

の

合

計
⑵

をフィルタリングの場合と同じ重みでリサンプリングすは固定ラグ平滑化
上述したとおり時点 t までの観測値 を与えらる こ と に よ っ て 得 ら れ る。 な お，
す る。

を得るには，mと時点 t －1における の一期先予測分布
れから粒子れた上で状態変数

を平滑化する場合，と仮定

固定点平滑化では，をから抜き出した

リサンプリングしてを生成

から生成された粒子である。北する。有限個の粒子（m 個）を用いて計算するた川（1996）によれば，この方法で固定区間平滑化

めに，t が大きな値の場合は何度もリサンプリン

4を実行できる。

グを繰り返すことによって，少数の粒子が繰り返4し出現してしまい，分布の形が損なわれることに北川（1996）の図６では、t が大きくなるにつれて固定
点平滑化分布が素早く固定区間平滑化分布－34－へ収束する様子が示されている。

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田洋志）なる。また，状態変数が k 次元ベクトルであれば，について平滑化のために k × t × m 個の粒子を保存しておについて以下のステップを実行すかなければならない。

を生成する。

る。

そこで，リサンプリングの対象となる粒子をであるからを生成する。（1）へ変更するととに基づいて生成したを用いるとともに，も に リ サ ン プ リ ン グ 回 数 を 減 ら す （t が小さという関係を利用して

い 場 合 は 増 加 ） こ と に よ っ て生成する。こうして一期先予測を行なう。

を生成する方法が直 前 で 計 算 し た

提案された。北川（2005）によれば，この方法でを 尤 度の式へ代入して，

固 定 ラ グ 平 滑 化 が 実 行 可 能 で あ る。 た だ し，を 生 成 す る。 こ のであり，L は固定した値とする。

からを計算するにはをを L 回リサンプリングすを 利 用 し て，をリサンプリングする。

こうしてフィルタリングを行なう。

る必要はあるが，t が大きくなった場合でもリサ尤度を 利 用 し て， 直 前 で 計 算 し た

ンプリングの回数が L 回に限定されるので退化と過去の平滑化の結果が進みにくい。これがこの方法の利点である。樋の合計である口編（2011）によれば，L は10～20に設定して平をリサンプリングすることに滑化を実行す る場合が多い。

5

この方法でもよ っ て 平 滑 化 が 実 行 さ れ て，から逐次的に粒子を生成していくことによってが得られる。

は得られる。そして，

はとを観測方程式（2）へ代2.6

入することによって得られる。

自己組織化状態空間モデル

パラメータ も状態変数として構築された状態平滑化する際に，固定点平滑化法では時点と同空間モデルは自己組織化状態空間モデルと呼ばれじ大きさの t 次のベクトルを生成していた。t ＞

る。このモデルを用いて未知のパラメータ と状

L+1の場合，固定ラグ平滑化では，t よりも小さ

態変数を同時に推定する。説明するに当たっては，い次元の L+1次のベクトルの生成で済むので，平

観測方程式の未知パラメータと遷移方程式の未知滑化のために保存しておかなければならない粒子

パラメータをそれぞれ と と区別して記述する。

の個数は k ×（L+1）× m 個となり，固定点平滑

つまり，

である。

化で実行するよりも少なくて済むため，計算上の遷移方程式

負荷が少ない。

観測方程式

（7）（8）

前述したとおり，自己組織化状態空間モデルでは，
2.5

モンテカルロフィルタのアルゴリズム

状態変数

モンテカルロフィルタにおける一期先予測とに加えて未知パラ

メータ を新たな状態変数と考え， を以下のとフィルタリング，平滑化（固定区間平滑化）を続おりに定義する。

けて行なうアルゴリズムは以下のように記述する



6

ことができる。

（9）を用いて記述した遷移方程式と観測方程式は以下のとおりである。

5

北川（1996）では L を20程度、最大でも50以下にするのが望ましいと記述されている。
6

ここでは確率分布から xt をサンプリングできると仮定しているが、サンプリングが難しい場合は提案分
布を用いて xt をサンプリングする。その場合は手遷移方程式

 （7’）

観測方程式

（8’）－35－

順が少し変更される。

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田

今までは未知パラメータ を一定値であると仮定洋志）る。 の確率密度関数は以下のとおりである。

していたが，そのままでは推定上の問題があるので， も時変すると仮定する。7

つまり， は以下のとおりに記述される。

（13）（10）と（13）から， の確率密度関数として次式を



（9’）

得る。

この場合，（7’）と（8’）を以下のように書き直すことができる。
において γ=1とおけば，パラメータ

遷移方程式

が1の指数分布となる。そこで， の確率分布と

観測方程式

（7”）

してパラメータが1の指数分布を仮定する。

（8”）



（14）

本節の2.2節から2.5節までに説明した方法を（10）と（14）から， の確率密度関数として次式用いて，（7”）と（8”）から構成される状態空間モデを得る。9ルにおける状態変数３を推定する。

SCD モデル

したがって，対数尤度は以下のとおりである。

以下では ACD モデルを説明した後に，SCD モ
デルを説明する。8

ただし， は（12a）と（12b）で定めたとおりである。
なお，林 ･ 佐藤（2016）によれば，ACD モデ

3.1

ACD モデルルの条件付き強度関数は以下のとおりである。

2つ事象が発生する時点 との間隔はと記述され，デュレーションと呼ばれる。
ACD モデルでは，観測される変数 が の条件付

ただし，はベースライン・ハザード関数であり，き期待値 （状態変数）と誤差項 を用いて以下の確率密度関数と生存関数の比である。

のように定式化される。


3.2

（10）


，あるいは を対数変換した

SCD モデル（11）

SCD モデルでは，観測される変数 が の条件の定式付き期待値 （状態変数）と誤差項 を用いて以化は様々なものを提案することが可能である。期下のように定式化される。

待デュレーション を定式化する基本モデルとして と の一次のラグ付き変数である

と



を説



（15）
（16）

明変数とする一次式を仮定する。



（12a）

このモデルではデュレーション の条件付き期待（12b）

値 が自らの過去の値

ここで，誤差項 の確率分布としてパラメータが γ で平均が１のガンマ分布
7を仮定すと誤差項 に依存している。
9であれば、であるから、畳み込みによって の確率密度関数は以下のとおりになる。

矢野（2014）を参照されたい。

ここで β=1を代入すると次式が得られる。

8

SCD モデルについては砂田（2018）においても説明さ
れている。

－36－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田期待デュレーション の定式化は自由である。

洋志）

変数 に加えて未知パラメータも状態変数として

また， の確率分布も0以上の値だけを実現する推定される。未知パラメータも推定の対確率分布という条件は必要であるものの，細かい象とするに当たり，条件はない。つまり， の定式化と の確率分布

ダムウォークに従うと仮定する。

は以下のようなランの特定化によって様々な SCD モデルを定義できる。
ポアソン過程では が独立かつ同一の指数分布に従っている。その結果，計数過程

分散 については非負制約が必要となるので対数がポア変換したものがランダムウォークに従うと仮定す

ソン過程に従っている。一方，SCD モデルでは，る。

が独立かつ同一の指数分布に従っていたとしても， が変動するので， の確率分布は変動する。
前述したとおり， ，あるいは を対数変換した

なお，誤差項は平均が０の正規分布に徒う。

の定式化は様々なものを提案することが可能である。本稿では単純なモデルとして，期新たな状態変数ベクトルを用いた遷移方程式は以下のとおりである。

待デュレーション を対数変換した の一階の自己回帰モデルを仮定する。

（18）

（17a）



（17b）

そして，新たな状態変数ベクトルを用いた観測方

このように SCD モデルでは誤差項を含むように

程式は以下のとおりである。

が定式化されるので，デュレーションに加えてボラティリティーにも確率的な変動が含まれるよ
うに ACD モデルが拡張されている。10
Bauwens=Veredas（2004）によれば，誤差項がワイブル分布に従う SCD モデルのハザード関



（19）

数（= 強度関数）は以下のとおりである。11ⅰ

一期先予測の方法（18）のの代わりに一期前のフィルタリング

した値

を代入して

の
3.3 自己組織化状態空間モデルを用いた推定の代わりに一期前のフィルタリングした値
を代入して

方法

を得る。次に，（18）を得る。そして，（18）のの代わりに一期前のフィルタリングした値

未知パラメータも状態変数として構築された状を代入して態空間モデルである自己組織化状態空間モデルを得る。さらに（18）,のの代わりに一期前のフィルタリングした値（Self-organizing state-space model）では，状態

の

を代入して

を得る。

10

ACD モデルでは を定式化する際に誤差項を含めない。
たとえば、と定式化する。

ⅱ

11

Bauwens=Veredas（2004） の389ペ ー ジ の（15）を 参 照されたい。ガンマ分布に従う場合のハザード関数は解析的に導出されないと記述されている。

－37－のフィルタリングの方法
尤度を重みとして利用して一期先予測の粒子からリサンプリングすることでモンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田フィルタリングした粒子
ⅲ

出来高
1
3
2
3
2

取引時刻
9：05：10
9：05：119：05：13… …

取引時刻
9：05：10
9：05：10
9：05：119：05：119：05：13を得る。

４のフィルタリングの方法

4.1

取引データ集計後
取引価格
1000
1002
1003

出来高
4
5
2
… …

取引データ集計前
取引価格
1000
1000
1002
1002
1003

取引データの集計例

… …

表1

洋志）

データとその加工方法について分析するデータについて

尤度を重みとして利用して一期先予測

本稿では東京商品取引所で取引されている金先の粒子からリサンプリングすることで物（ 現 物 受 渡 用 ） の テ ィ ッ ク デ ー タ を 用 い てフィルタリングした粒子
尤度を得る。同様に，と一期先予測の粒子てフィルタリングした粒子に，尤度

SCD モデルを推定した。試験的に2016年4月11を用い日を選んでいる。金先物も含めて，日本の商品先を得る。さら

物市場では一般に期先物の出来高が最も多い。そと一期先予測の粒子

こで，今回の分析でも4月11日における期先物でを得る。

ある2017年2月限のデータを利用した。金先物市を利用して粒子場における取引時間は午前9時～15時15分が日中
ⅳの平滑化の方法（固定点平滑化）
尤度の取引時間である。9時に板合せ取引が行なわれを重みとして利用して，一期前のて9時1分に始値が付く。その後，ザラバ取引が平滑化で生成した粒子

15時15分まで行なわれる。12とに時間外取引も行われているが，本稿の分析ではの合計からリサンプリングすることで平滑化した粒子を得る。

これ以外の時間帯

9時1分から15時15分までの日中のデータだけを利用してパラメータの推定を行った。分析対象とⅴ

の平滑化の方法（固定点平滑化）

尤度を重みとして利用して，一期前の平滑化で生成した粒子で取引される一つ一つの取引に関するデータを集の合計をリサンプリングすることで平滑めたものである。したがって，ティックデータに尤度を得る。同様に，は，取引時刻（時分秒），取引価格（円），取引数と一期前の平滑化で生成した粒子

量（枚）が含まれている。東京商品取引所から購との合計を用いて平滑入したティックデータには，日付（年月日），取を得る。さらに，

引種別（先物かオプションか），商品番号，限月（年と一期前の平滑化で生成した粒子月），約定時刻（時分秒），約定価格（円），出来

化した粒子
尤度

と
い

本稿で分析するティックデータとは，1日の中

と

化した粒子

用

する取引は4，017個であった。

て

平

滑

化

の合計を
し

た

を得る。

粒

高（枚）が入力されている。

子

ティックデータを用いて分析するので，1秒間に複数の取引が執行される場合もある。この場合
12

現在は取引システムが変更されて、8時45分から15時15分までが日中の取引時間であり、寄り付きと大引けで板合せ取引を行なう。

－38－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田図１

デュレーション洋志）

を平滑した結果（縦軸は取引間隔，横軸は経過時間）は表1に示すとおり，同じ時点に複数の取引情報

具体的には，誤差項に指数分布を仮定した SCDが表示される。取引間隔を計算する場合には同じ

モデル，つまり，（15），（17a，b）と（14）式に基時点（秒）で執行された取引を分けることはできづいたモデルのパラメータを自己組織化状態空間ないので，1つのデータと考える。こうして，取モデル（18），（19）として推定した。推定に当たっ引が執行された時点（秒）ごとに出来高（取引枚ては，各取引時点で5，000個の粒子を発生させた。

数）と約定価格が2，088個得られた。

未知パラメータを推定する際の初期値として SCD モデルをベイズ推定した結果（表3）4.2

データの加工方法についてを利用している。また，未知パラメータはランダ

全ての取引が含まれているティックデータから，

ム ウ ォ ー ク に 従 う と 仮 定 し た が， 誤 差 項

時間間隔データを作成する。そのデータにSCD モ
デルを直接適用して，パラメータを推定するのでの標準偏差は杉浦 = 中妻（2016）
に従い，0.05と仮定した。

はなく，時間間隔データを元に計算した平滑化曲

推定を行なう前に取引間隔（デュレーション）

線から得られた値でデータを割った後，そのデーのデータを平滑化によって加工しなくてはいけな

タを用いてSCD モデルのパラメータを推定する。

い。平滑化の方法は幾つかあるが，森保（2006）

平滑化の方法は様々である。局所多項式回帰にに従ってスーパースムーザを利用して，取引間隔

おいてバンド幅を可変とした LOESS（LOcallyのデータを平滑化した。原データを平滑化した値

weighted RegrESSion） が 考 案 さ れ た。LOESS

で割ることによって確定的な日中の傾向を予め除をさらに発展させたものにスーパースムーザがあ去しておくのである。14る。本稿ではスーパースムーザを用いてデータを結果が図1である。横軸は取引開始時刻近く（913

平滑化する。

取引間隔を平滑化した時）からの経過時間（単位は秒）である。午前中

５は最初の1時間半，次の1時間，午後に山があり，推定結果

大きく3つの山があることが読み取れる。

本稿では試験的に2016年4月11日に東京商品取

平滑化した値で取引間隔 を割って調整した値

引所で取引された金先物市場（現物受渡用）の
14

Engle=Russell（1998）の1137ページに期待デュレーショ
ンは確定的（deterministic）な成分と確率的（stochastic）
な成分に分けられると記述されている。

ティックデータを用いて SCD モデルを推定した。
13

データの平滑化については竹澤（2007）に詳しい。

－39－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田

図２

図３

表２

洋志）

の推移（縦軸は平均，横軸は取引された順番），，の推移（縦軸は平均，横軸は取引された順番）

SCDモデルを自己組織化状態空間モデルとして推定した結果
変数

全取引時点の平均

全取引時点の標準偏差－0.001833

0.108811

0.862557

0.098860

0.026804

0.003024－40－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田

表3
パラメータ

洋志）

SCDモデルをベイズ推定した結果
事後平均

事後標準偏差事後自己相関0.004955

0.005017

0.220067

0.855059

0.022434

0.409312

0.027147

0.002642

0.712333
出典：砂田（2018）を改めてデータと定めて SCD モデルを当ては均が－0.0018，標準偏差が0.108811であった。各める。自己組織化状態空間モデルを用いるので，取引時点において生成した5，000個の粒子の標本状態変数だけでなく，パラメータを推定した結果

平均を標本標準偏差で割った値の絶対値は全てのも取引時点ごとに計算される。図2に状態変数の取引時点で1.96以下であった。したがって，有意系列に0と異なるとは判断できないであろう。なお，をフィルタリングによって推定した結果 を 示 す。

15

各 取 引 時 点 で5，000個 の 粒 子を生成してその平均を各取引時点にであるから，定数項 が負値になっても問題は無い。回帰係数 については，の全おける状態変数の推定値と考える。図2の縦軸は取 引 時 点 を 通 じ た 平 均 が0.8626， 標 準 偏 差 が各取引時点で状態変数を5，000回生成した結果を0.09886であるから大きな変化がないと考えられ平均した値る。各取引時点において生成した5，000個の粒子，横軸は取引された順番である。

図2を見ると，の系列は0の周囲を変動していることが分かる。

の標本平均を標本標準偏差で割った値は全ての取であるから の推定値

引時点で1.96以上であった。誤差項の分散 につは多くの取引時点において1に近い値であると考

いては，

えられる。

0.0268，標準偏差が0.0030であるから大きな変化

次に，図3に3つのパラメータの全取引時点を通じた平均がをフィがないと考えられる。各取引時点において生成しルタリングによって推定した結果を示す。3つのた5，000個の粒子の標本平均を標本標準偏差でパラメータのそれぞれについて，各取引時点で割った値は全ての取引時点で1.96以上であった。

5，000個の粒子して，その平均（ ，を生成，

また，DIC は4453.203，）を各取引時点にを用いて計算した対数尤度は－1935.28であった。

おけるパラメータの推定値と考えている。図3の参考のため，ベイズ推定した結果を表3に示し

縦軸は各取引時点で3つのパラメータをそれぞておく。16

れ5，000回生成した結果を平均した値，横軸は取であることを棄却できない。しかし，回帰係数引された順番である。3つのパラメータのそれぞは0.855であり有意に0と異なると判定された。

れについて2，088個の推定値があり，その推移を採択率は約97.7% であった。

示したところ，2，088個の取引時点を通じて3つ

56.7% であった。また，DIC は4472.84，状態変数ともほとんど変化していないように見える。

の系列

表2には，3つのパラメータの2，088個の採択率は約の事後平均を用いて計算した対数尤

度は－2117.84であった。以上から ， と の推定の推定値（= 平均）について全取引時間を通じて

値はおおよそ似た結果となった。

計算した標本平均と標本標準偏差を示す。定数項については，

定数項 は0.00496であり，有意に0６の全取引時点を通じた平結論

本稿の目的はモンテカルロフィルタを紹介する

15

2節と3節において平滑化について説明したが、平滑化による推定までは行なわない。

16

－41－

砂田（2018）の推定結果を再掲したものである。

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田洋志）

ことと，モンテカルロフィルタを用いて SCD モ［6］林 高樹，2012，「高頻度データの分析⑵：不デルを推定する方法を紹介することである。第2等間隔データ分析」，刈屋武昭他編，『経済時節では，自己組織化状態空間モデルを含めてモン系列ハンドブック』の第6章4節に所収，朝テカルロフィルタを紹介した後，その推定方法を倉書店 .

紹介した。第3節では，SCD モデルと自己組織［7］樋口知之編，2011年，『データ同化入門』，朝化状態空間モデルによる推定方法を説明した。第倉書店．

4節では，分析対象の取引間隔（デュレーション）［8］樋口知之，2011年，『予測に生かす統計モデ

データについて説明した後，日中効果を取り除く

リングの基本』，講談社．

のに用いる平滑化について説明した。平滑化の手［9］Kalman，R.E.，1960，“A new approach to

法はいろいろあるが，本稿では森保（2006）と同

linear filtering and prediction

じスーパースムーザを利用して平滑化したデータ

problems，”Trans. ASME-J，Basic Eng.（ser.D），を利用して推定した。第5節では，東京商品取引82，pp.35-45.

所の金先物市場のティックデータを対象にして［10］Kitagawa， G.， 1993， “A monte carlo

SCD モデルのパラメータをモンテカルロフィル

filtering and smoothing method for non-

タで推定した結果を示した。SCD モデルをベイ

Gaussian nonlinear state space

ズ推定した結果と比較したところ，類似の結果がmodels，”Technical report， Institute of得られた。

Statistical Mathematics researchmemorandum，No.462.

参考文献［11］Kitagawa，G.，1996，“Monte carlo filter and［1］Bauwens，L.，and Veredas，D.，2004，“Thesmoother for non-Gaussian nonlinear state

stochastic conditional duration model: aspace models，”Journal of Computational andlatent variable model for the analysis ofGraphical Statistics，5，pp.1-25.

financial durations，”Journal of Econometrics，［12］北川源四郎，“ モンテカルロ ･ フィルタおよ119，pp.381-412.

び平滑化について ”，1996，『統計数理』，第［2］E ngle， R.F.， and Russel， J.R.， 1997，

44巻1号，pp.31-48．

“Forecasting the frequency of changes in［13］K itagawa， G.， 1998， “A self-organizing

quoted foreign exchange price with thestate space models，”Journal of American

ACD model，”Journal of Empirical Finance，

Statistical Association，93，pp.1203-1205.

4 ，pp.187 - 212.

［14］北川源四郎，2005年，『時系列解析入門』，岩［3］E ngle， R.F. and Russel， J.R.， 1998，“Autoregressive conditional duration: a new

波書店．
［15］Men，Z.，and Kolkiewicz，A. W.，Wirjanto，model for irregularly spaced transaction

T . S . ， 2015 ， “ B a y e s i a n a n a l y s i s o fdata，”Econometrica，66，pp.1127-1162.

asymmetric stochastic conditional duration［4］Gordon，N.，Salmond，D. and Smith，A.，

model，”Journal of Forecasting，34，pp.36 -1993，“Novel approach to nonlinear /nonGaussian Bayesian state estimation，”IEEE

56.
［16］Men，Z.，and Kolkiewicz，A. W.，Wirjanto，Proceedings-F，140，pp.107-113.

T . S . ， 2016 ， “ B a y e s i a n i n f e r e n c e o f［5］林 高樹・佐藤彰洋，2016，『金融市場の高頻asymmetric stochastic conditional duration

度データ分析』，朝倉書店 .

models，”Journal of Statistical Computation－42－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田and Simulation，86，pp.1295 - 1319.

洋志）

金利制約下における日本経済－流動性制約家［17］森保洋，「金先物市場の日中価格変動と取引計を含むニューケインジアン DSGE モデル」，

時間間隔」，2006，『経営と経済』，長崎大学『世界同時不況と景気循環分析』の第9章に経済学会，第85巻3･4号，pp.408-430.

所収，浅子 ･ 飯塚 ･ 宮川編集，東京大学出版［18］P acurar， M.， 2008，“Autoregressive

会.

conditional duration models in finance : A

2014，“ 粒子フィルタの基礎と応用：［27］矢野浩一，

survey of the theoretical and empirical

フィルタ・平滑化・パラメータ推定 ”，『日literature，”Journal of Economic Survey，22，

本統計学会誌』，第44巻，pp.189-216.

pp.711-751.
［19］佐 藤整尚 ･ 高橋明彦，2005，“ モンテカルロ
フィルタを用いた金利モデルの推定 ”『計算統計Ⅱ－マルコフ連鎖モンテカルロ法の基礎』の補論 B に所収，岩波書店 .
［20］砂 田洋志，2018，“SCD モデルの理論と実証―東京商品取引所における金先物取引への応用―”，山形大学紀要（社会科学），第48巻第2号，pp.41-63.
［21］S trickland， C.M.， Forbes， C.S.， andMartin，G.M.，2006，“Bayesian analysis ofstochastic conditional duration
model，”Computational Statistics & Data
Analysis，50，pp.2247-2267.
［22］杉 浦航 ･ 中妻照雄，2016，「ティックデータを用いた株式市場における約定予測」，日本
金融・証券計量・工学学会編『ファイナンス
における数値計算手法の新展開』の第4章に所収，朝倉書店 .
［23］高部勲，2010，“ 非線形状態空間モデル及びモンテカルロ ･ フィルタによる NAIRU の推定 ”，『統計研究彙報』，第67号，pp.55-88．
［24］竹 澤邦夫，2007，『みんなのためのノンパラ
メトリック回帰（上）第3版』，吉岡書店 .
［25］矢 野浩一・佐藤整尚，2007，“ 初期分布探索付き自己組織化状態空間モデルによる金融時
系列解析の最前線：t 分布付き確率的ボラ
ティリティ変動モデルへの応用 ”，『FSA リ
サーチレビュー 2006』，金融庁金融研究セン
ター，pp1-24.
［26］矢 野浩一 ･ 飯田泰之 ･ 和合肇，2011，
「ゼロ
－43－

モンテカルロフィルタの理論と応用――SCD モデルへの応用――（砂田洋志）

Theory and Application of Monte Carlo Filter―― Application to Stochastic Conditional Duration model ――
Hiroshi SUNADA（Management, Social Systems Course）
The monte carlo filter is a very useful estimation method, because we can use it even if the error
term does not follow the normal distribution and /or model is non-linear. In this paper we introduce
monte carlo filter and show how to apply it to financial model. We apply it to stochastic conditional
duration model. We use tick data of golden futures on Tokyo Commodity Exchange. We estimate theparameters and state variables by monte carlo filter. The result is almost the same as that by Bayesian
method.

－44－