R言語 パッケージ {sde} をつかって確率微分方程式モデル

  • 7
    いいね
  • 0
    コメント
この記事は最終更新日から1年以上が経過しています。

R Package sde をつかってみる

【参考文献】

Iacus, Stefano M.(著) "Simulation and Inference for Stochastic Differential Equations"

Book.jpg

刊行元 Springer社 公式ウェブサイトへのリンク

【以下、 sde のhelp()ページマニュアルを参考に使ってみる】

パッケージの読込み

require("sde")

ブラウン運動系列を発生させてみる

【関数】 BM() ※ Brownian motion関数

※ブラウン運動 変数の発生シミュレーション関数

ブラウン運動の引数パラメータ群と今回の設定

※関数を実行するたびに、発生されるブラウン(確率過程)時系列値は変化する

  • 初期値 $x$ : 0
  • 初期時刻 $t_0$ : 0
  • 最終時刻 $T$ : 1
  • 各時系列変数の数 $N$ : 1,000

※(時系列変数の数 = グラフ横軸の刻み幅 )

  • ( $T$ - $t_0$ )/ $N$

実行(1回目)

BMpath <- BM(x = 0, t0 = 0, T = 1, N = 1000) 

plot(BMpath)

graph.1.png

実行(2回目)

BMpath <- BM(x = 0, t0 = 0, T = 1, N = 1000) 

plot(BMpath)

graph.2.png

実行(3回目)

BMpath <- BM(x = 0, t0 = 0, T = 1, N = 1000) 

plot(BMpath)

graph.3.png

株価過程のダミーとして使われる幾何ブラウン運動系列を発生させてみる

【関数】 GBM() ※ Geometric Brownian motion関数

※幾何ブラウン運動 変数の発生シミュレーション関数

【公式】


dx(t)=\mu x(t)dt+\sigma x(t)dW_t

【GBM()関数の返り値】


X(t)=x \mathrm{e}^\frac {r- \sigma^2} {2}*t+\sigma W_t

【解く式】


X(t)=x \mathrm{e}^\frac {0 - 1^2} {2}*t+1 W_t

※関数を実行するたびに、(確率過程)のシミュレーション結果(時系列値の推移)は変化する

-(※返り値) 観測不能な時系列変数 (an invisible ts object): $x$
- 最終時刻 $T$ : 1
- 初期値 $x$ : 1
- 金利 $r$ : 0
- ボラティリティ $σ$ : 1
- 各時系列変数の数 $N$ : 1,000

※(時系列変数の数 = グラフ横軸の刻み幅 )

  • ( $T$ - $t_0$ )/ $N$

実行(1回目)

GBMpath <- GBM(x=1, r=0, sigma=1, T=1, N=100)

plot(GBMpath)

graph.4.png

実行(2回目)

GBMpath <- GBM(x=1, r=0, sigma=1, T=1, N=100)

plot(GBMpath)

graph.5.png

実行(3回目)

GBMpath <- GBM(x=1, r=0, sigma=1, T=1, N=100)

plot(GBMpath)

graph.6.png

helpページのサンプル例に沿って確率微分方程式を解いてみる

(解の時系列推移の概形をグラフで確認)

【関数】 sde.sim()

(stochastic differential equation simulation関数)

※返り値: N+1個(時点)の解のベクトル(ts型オブジェクト)

(注)上記 N は、関数の引数に与えて設定する(以下参照)

【sde.sim()関数の書式】

sde.sim(t0 = 0, T = 1, X0 = 1, N = 100, 100, drift, sigma, drift.x, sigma.x, drift.xx, sigma.xx, drift.t, method = c("euler", "milstein", "KPS", "milstein2", "cdist","ozaki","shoji","EA"), alpha = 0.5, eta = 0.5, pred.corr = T, rcdist = NULL, theta = NULL, model = c("CIR", "VAS", "OU", "BS"), k1, k2, phi, max.psi = 1000, rh, A, M=1)

  • 初期時刻 : $t_0$
  • 最終時刻 : $T$
  • 初期値 : $x_0$
  • シミュレーション回数(確率変数の個数): $N$
  • 出力する時系列推移(trajectory)の数: $M$
  • シミュレーションのTime Step: $δ$
  • ドリフト項の係数(確定項): drift
  • 拡散項係数(確率項):sigma
  • ドリフト項係数の偏微分値:drift.x
  • 拡散項係数の偏微分値: sigma.x
  • ドリフト項係数の2次偏微分値: drift.xx
  • 拡散項係数の2次偏微分値: sigma.xx
  • ドリフト時間項係数の偏微分値: drift.t
  • シミュレーション演算法: method

(以下略)

オルンシュタイン=ウーレンベック過程(「平均回帰過程」) (Ornstein-Uhlenbeck process)

【公式】


dr(t)=-\theta(r(t)-\mu)dt+\sigma dW_t

【解く式】


dX(t)=-5(X(t)-0)dt+3.5 dW_t

  • $θ$ : 5
  • $r$ : X
  • $μ$ : 0
  • $σ$ : 3.5

(Wikipedia)オルンシュタイン=ウーレンベック過程

「オルンシュタイン=ウーレンベック過程は、離散時間AR(1)過程の連続時間バージョンであると言える」

set.seed(123)

d <- expression(-5 * x)   # (注意)小文字の"x"
s <- expression(3.5) 

X <- sde.sim(X0=10,drift=d, sigma=s)  # (注意)大文字の"X"

plot(X,main="Ornstein-Uhlenbeck")   # (注意)大文字の"X"

graph.7.png

Multiple trajectories of the O-U process

※ 引数 M に3を与えて、微分方程式の解として、3本の時系列(過程)を得ている

set.seed(123)

X <- sde.sim(X0=10,drift=d, sigma=s, M=3) # (注意)大文字の"X"

plot(X,main="Multiple trajectories of O-U")  # (注意)大文字の"X"

graph.8.png

※(ケース2)引数 M に5を与えた場合

set.seed(123)

X <- sde.sim(X0=10,drift=d, sigma=s, M=5) 

plot(X,main="Multiple trajectories of O-U\n (M=5 )") 

addition.graph.png

コックス・インガーソル・ロス・モデル (CIRモデル: Cox-Ingersoll-Ross process)

【公式】


dr_t=a(b-r_t)dt+\sigma \sqrt{r_t} dW_t

【解く式】


dX_t=3(2-X_t)dt+2 \sqrt{X_t} dW_t

  • $r$ : X
  • $a$ : 3
  • $b$ : 2
  • $σ$ : 2

(Wikipedia)コックス・インガーソル・ロス・モデル

「数理ファイナンスにおいて利子率の時間的変動を記述する数理モデルの一つ」

set.seed(123)

d <- expression( 6-3*x ) 
s <- expression( 2*sqrt(x) ) 

X <- sde.sim(X0=10,drift=d, sigma=s)   # (注意)大文字の"X"

plot(X,main="Cox-Ingersoll-Ross")  # (注意)大文字の"X"

graph.9.png

Cox-Ingersoll-Ross using the conditional distribution "rcCIR"

set.seed(123)

X <- sde.sim(X0=10, theta=c(6, 3, 2), rcdist=rcCIR, method="cdist") 

plot(X, main="Cox-Ingersoll-Ross")

graph.10.png

Exact simulation ※同義語 "perfect simulation" , "perfect simulation", CFTP(coupling from the past)

(参考)朱鷲の杜Wiki

Markovモデルの定常状態でのサンプリングする方法.
MCMCによるサンプリングは近似だが,パーフェクトサンプリングは真の同時分布に従うサンプルを得るように設計されている.でも,計算はずっと大変.
過去からのカップリングが代表的な手法.

set.seed(123)

d <- expression(sin(x))
d.x <- expression(cos(x)) 

A <- function(x) 1-cos(x)

X <- sde.sim(method="EA", delta=1/20, X0=0, N=500, drift=d, drift.x = d.x, A=A)

plot(X, main="Periodic drift")

graph.11.png

(参考ウェブサイト)

(参考 YouTube 動画(※ 英語))

上記の動画が第1回目の10本以上の動画シリーズ(英語)です

動画投稿者は、Quant Education さんです。

統計数値計算・時系列解析など、R やその他のツールを用いた実装方法の使い方チュートリアル画面を提供してくれています

YouTube Quant Education 投稿動画リストへのリンク )

I have worked as a Senior Quant for a big bank and studied at London School of Economics, Oxford University and Cambridge University.
The purpose of this channel is to widely disseminate the knowledge of tools and techniques used in Mathematical/Quantitative Finance and Quantitative Business Analytics. Please also visit QuantEdu.com for more courses. All the material here is taken from textbooks and is standard approach, my judgement and feeling is not part of the lectures.
I also do consulting and onsite training for R, Probability, Time Series and Monte Carlo simulation. I am based in Pennsylvania, USA. Contact (training/partnership) :: quantedu.com@gmail.com or leave me a message on this channel.