library(sde) (R)

Rのパッケージ，library(sde)の紹介です．

このパッケージでは確率過程のシミュレーションが簡単にできます．

たとえばブラウン運動を発生させたければ

> BM
function (x = 0, t0 = 0, T = 1, N = 100) 

という関数を使うことで，
初期値$x$，初期時刻$t_0$，最終時刻$T$，シミュレーション刻み幅$(T-t_0)/N$
のブラウン運動のパスが生成できます．

この関数を使って

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

とすれば，下のようなブラウン運動が生成できます．

このほかにも，ブラウン橋

BBridge(x=0, y=0, t0=0, T=1, N=100)

や確率微分方程式の解

sde.sim(t0 = 0, T = 1, X0 = 1, N = 100, delta, 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)

のシミュレーションなども簡単にでき，便利です．

最後に$B_t$を0を出発する標準ブラウン運動として
\[
\mathrm{E} \bigg[B_t\int_0^t B_s {\rm d}s \bigg],
\]
の値をモンテカルロ法で計算したいと思います．

ちなみに解析解は，まず伊藤の公式から
\[
\int_0^t B_s {\rm d}s=tB_t-\int_0^t s {\rm d}B_s,
\]
なので
\begin{eqnarray*}
\mathrm{E} \bigg[B_t\int_0^t B_s {\rm d}s \bigg]&=&t\mathrm{E}[B_t^2]-E[B_t\int_0^t s {\rm d}B_s]\\
&=&t^2-\mathrm{E}\langle B_t, \int_0^t s {\rm d}B_s\rangle\\
&=&t^2-\mathrm{E}\int_0^t s {\rm d}s\\
&=&\frac{t^2}{2}.
\end{eqnarray*}

Rのコードは下の通り．

#E[B_t \int_0^t B_s ds]の計算（モンテカルロ法）
library(sde)
time_max <- 10 #[0,time_max]までBMをシミュレーション
sim_num <- 1000 #BMの各時刻での発生回数
crossvar_mean <- rep(NaN, time_max)
crossvar <- rep(NaN, sim_num)
for(i in 1:time_max){
	for(j in 1:sim_num){
		#BMのpathを生成
		BMpath <- BM(T=i, N=100*i)

		#\int_0^t B_s ds (リーマン和)
		int_Bs_ds <- sum(BMpath[1:(100*i)]/100) 

		crossvar[j] <- BMpath[100*i+1]*int_Bs_ds
	}
	crossvar_mean[i] <- mean(crossvar)
}

plot(crossvar_mean)

結果は以下の通り．

解析解と一致していることが分かります．