# ぱらぱらめくる『Simulation and Inference for Stochastic Differential Equations: With R Examples』

Simulation and Inference for Stochastic Differential Equations: With R Examples (Springer Series in Statistics)

# ぱらぱらめくる『Simulation and Inference for Stochastic Differential Equations: With R Examples』Rソース

• Rにはsdeというパッケージがあって、それのモト本が、この本らしい
• ちょっと重い処理も含まれるが、examplesをひたすら抜き出すと以下の通り
library(sde)
plot(BM())
plot(BBridge())
plot(GBM())

tau0 <- 0.6
k0 <- ceiling(1000*tau0)
set.seed(123)
X1 <- sde.sim(X0=1, N=2*k0, t0=0, T=tau0, model="CIR",
theta=c(6,2,1))
X2 <- sde.sim(X0=X1[2*k0+1], N=2*(1000-k0), t0=tau0,
T=1, model="CIR", theta=c(6,2,3))
Y <- ts(c(X1,X2[-1]), start=0, deltat=deltat(X1))
X <- window(Y,deltat=0.01)
DELTA <- deltat(X)
n <- length(X)
mu <- function(x) 6-2*x
sigma <- function(x) sqrt(x)
cp <- cpoint(X,mu,sigma)
cp
plot(X)
abline(v=tau0,lty=3)
abline(v=cp$tau0,col="red") # nonparametric estimation cpoint(X) d <- expression((3-x)) s <- expression(1.2*sqrt(x)) par(mar=c(3,3,1,1)) par(mfrow=c(2,1)) set.seed(123) X <- DBridge(x=1.7,y=0.5, delta=0.01, drift=d, sigma=s) plot(X) X <- DBridge(x=1,y=5, delta=0.01, drift=d, sigma=s) plot(X) ## Not run: d1 <- function(t,x,theta) theta[1]*(theta[2]-x) s1 <- function(t,x,theta) theta[3]*sqrt(x) from <- 0.08 x <- seq(0,0.2, length=100) sle10 <- NULL sle2 <- NULL sle5 <- NULL true <- NULL set.seed(123) for(to in x){ sle2 <- c(sle2, dcSim(from, to, 0.5, d1, s1, theta=c(2,0.02,0.15), M=50000,N=2)) sle5 <- c(sle5, dcSim(from, to, 0.5, d1, s1, theta=c(2,0.02,0.15), M=50000,N=5)) sle10 <- c(sle10, dcSim(from, to, 0.5, d1, s1, theta=c(2,0.02,0.15), M=50000,N=10)) true <- c(true, dcCIR(to, 0.5, from, c(2*0.02,2,0.15))) } par(mar=c(5,5,1,1)) plot(x, true, type="l", ylab="conditional density") lines(x, sle2, lty=4) lines(x, sle5, lty=2) lines(x, sle10, lty=3) legend(0.15,20, legend=c("exact","N=2", "N=5", "N=10"), lty=c(1,2,4,3)) ## End(Not run) data(DWJ) ret <- diff(DWJ)/DWJ[-length(DWJ)] par(mfrow=c(2,1)) par(mar=c(3,3,2,1)) plot(DWJ,main="Dow-Jones closings",ylab="",type="p") plot(ret,main="Dow-Jones returns",ylab="",type="p") cp <- cpoint(ret) cp abline(v=cp$tau0,lty=3)
cp <- cpoint(window(ret,end=cp$tau0)) cp abline(v=cp$tau0,lty=3)

set.seed(123)
d <- expression(-1*x)
s <- expression(2)
sde.sim(drift=d, sigma=s) -> X
S <- function(t, x, theta) sqrt(theta[2])
B <- function(t, x, theta) -theta[1]*x
true.loglik <- function(theta){
DELTA <- deltat(X)
lik <- 0
for(i in 2:length(X))
lik <- lik + dnorm(X[i], mean=X[i-1]*exp(-theta[1]*DELTA),
sd = sqrt((1-exp(-2*theta[1]*DELTA))*
theta[2]/(2*theta[1])),TRUE)
lik
}
xx <- seq(-3,3,length=100)
sapply(xx, function(x) true.loglik(c(x,4))) -> py
sapply(xx, function(x) EULERloglik(X,c(x,4),B,S)) -> pz
# true likelihood
plot(xx,py,type="l",xlab=expression(beta),ylab="log-likelihood")
lines(xx,pz, lty=2) # Euler

## Not run:
alpha <- 0.5
beta <- 0.2
sigma <- sqrt(0.05)
true <- c(alpha, beta, sigma)
names(true) <- c("alpha", "beta", "sigma")
x0 <- rsCIR(1,theta=true)
set.seed(123)
sde.sim(X0=x0,model="CIR",theta=true,N=500000,delta=0.001) -> X
X <- window(X, deltat=0.1)
DELTA = deltat(X)
n <- length(X)
X <- window(X, start=n*DELTA*0.5)
plot(X)
u <- function(x, y, theta, DELTA){
c.mean <- theta[1]/theta[2] +
(y-theta[1]/theta[2])*exp(-theta[2]*DELTA)
c.var <- ((y*theta[3]^2 *
(exp(-theta[2]*DELTA)-exp(-2*theta[2]*DELTA))/theta[2] +
theta[1]*theta[3]^2*
(1-exp(-2*theta[2]*DELTA))/(2*theta[2]^2)))
cbind(x-c.mean,y*(x-c.mean), c.var-(x-c.mean)^2,
y*(c.var-(x-c.mean)^2))
}
CIR.lik <- function(theta1,theta2,theta3) {
n <- length(X)
dt <- deltat(X)
-sum(dcCIR(x=X[2:n], Dt=dt, x0=X[1:(n-1)],
theta=c(theta1,theta2,theta3), log=TRUE))
}
fit <- mle(CIR.lik, start=list(theta1=.1, theta2=.1,theta3=.3),
method="L-BFGS-B",lower=c(0.001,0.001,0.001), upper=c(1,1,1))
# maximum likelihood estimates
coef(fit)
gmm(X,u, guess=as.numeric(coef(fit)), lower=c(0,0,0),
upper=c(1,1,1))
true
## End(Not run)

set.seed(123)
d <- expression(-1*x)
s <- expression(2)
sde.sim(drift=d, sigma=s) -> X
M0 <- function(t, x, theta) -theta[1]*x
M1 <- function(t, x, theta) -theta[1]
M2 <- function(t, x, theta) 0
M3 <- function(t, x, theta) 0
M4 <- function(t, x, theta) 0
M5 <- function(t, x, theta) 0
M6 <- function(t, x, theta) 0
mu <- list(M0, M1, M2, M3, M4, M5, M6)
F <- function(t, x, theta) x/sqrt(theta[2])
S <- function(t, x, theta) sqrt(theta[2])
true.loglik <- function(theta) {
DELTA <- deltat(X)
lik <- 0
for(i in 2:length(X))
lik <- lik + dnorm(X[i], mean=X[i-1]*exp(-theta[1]*DELTA),
sd = sqrt((1-exp(-2*theta[1]*DELTA))*theta[2]/
(2*theta[1])),TRUE)
lik
}
xx <- seq(-3,3,length=100)
sapply(xx, function(x) HPloglik(X,c(x,4),mu,F,S)) -> px
sapply(xx, function(x) true.loglik(c(x,4))) -> py
plot(xx,px,type="l",xlab=expression(beta),ylab="log-likelihood")
lines(xx,py, lty=3) # true

set.seed(123)
theta <- c(6,2,1)
X <- sde.sim(X0 = rsCIR(1, theta), model="CIR", theta=theta,
N=1000,delta=0.1)
b <- function(x)
theta[1]-theta[2]*x
sigma <- function(x)
theta[3]*sqrt(x)
minX <- min(X)
maxX <- max(X)
par(mfrow=c(3,1))
curve(b,minX,maxX)
lines(ksdrift(X),lty=3)
curve(sigma,minX, maxX)
lines(ksdiff(X),lty=3)
f <-function(x) dsCIR(x, theta)
curve(f,minX,maxX)
lines(ksdens(X),lty=3)

set.seed(123)
d <- expression(-1 * x)
s <- expression(1)
x0 <- rnorm(1,sd=sqrt(1/2))
sde.sim(X0=x0,drift=d, sigma=s,N=1000,delta=0.1) -> X
d <- expression(-theta * x)
linear.mart.ef(X, d, s, a1=expression(-x), lower=0, upper=Inf,
c.mean=expression(x*exp(-theta*0.1)),
c.var=expression((1-exp(-2*theta*0.1))/(2*theta)))

data(quotes)
plot(quotes)
d <- MOdist(quotes)
cl <- hclust( d )
groups <- cutree(cl, k=4)
cmd <- cmdscale(d)
plot( cmd, col=groups)
text( cmd, labels(d) , col=groups)
plot(quotes, col=groups)
plot(quotes, col=groups,ylim=range(quotes))

data(quotes)
plot(quotes)
data(quotes)
plot(quotes)

rcBS(n=1, Dt=0.1, x0=1, theta=c(2,1))

rcCIR(n=1, Dt=0.1, x0=1, theta=c(6,2,2))

rcOU(n=1, Dt=0.1, x0=1, theta=c(0,2,1))

rsCIR(n=1, theta=c(6,2,1))

rsOU(n=1, theta=c(0,2,1))

# Ornstein-Uhlenbeck process
set.seed(123)
d <- expression(-5 * x)
s <- expression(3.5)
sde.sim(X0=10,drift=d, sigma=s) -> X
plot(X,main="Ornstein-Uhlenbeck")
# Multiple trajectories of the O-U process
set.seed(123)
sde.sim(X0=10,drift=d, sigma=s, M=3) -> X
plot(X,main="Multiple trajectories of O-U")
# Cox-Ingersoll-Ross process
# dXt = (6-3*Xt)*dt + 2*sqrt(Xt)*dWt
set.seed(123)
d <- expression( 6-3*x )
s <- expression( 2*sqrt(x) )
sde.sim(X0=10,drift=d, sigma=s) -> X
plot(X,main="Cox-Ingersoll-Ross")
# Cox-Ingersoll-Ross using the conditional distribution "rcCIR"
set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), rcdist=rcCIR,
method="cdist") -> X
plot(X, main="Cox-Ingersoll-Ross")
set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), model="CIR") -> X
plot(X, main="Cox-Ingersoll-Ross")
# Exact simulation
set.seed(123)
d <- expression(sin(x))
d.x <- expression(cos(x))
A <- function(x) 1-cos(x)
sde.sim(method="EA", delta=1/20, X0=0, N=500,
drift=d, drift.x = d.x, A=A) -> X
plot(X, main="Periodic drift")

set.seed(123)
# true model generating data
dri <- expression(-(x-10))
dif <- expression(2*sqrt(x))
sde.sim(X0=10,drift=dri, sigma=dif,N=1000,delta=0.1) -> X
# we test the true model against two competing models
b <- function(x,theta) -theta[1]*(x-theta[2])
b.x <- function(x,theta) -theta[1]+0*x
s <- function(x,theta) theta[3]*sqrt(x)
s.x <- function(x,theta) theta[3]/(2*sqrt(x))
s.xx <- function(x,theta) -theta[3]/(4*x^1.5)
# AIC for the true model
sdeAIC(X, NULL, b, s, b.x, s.x, s.xx, guess=c(1,1,1),
lower=rep(1e-3,3), method="L-BFGS-B")
s <- function(x,theta) sqrt(theta[3]*+theta[4]*x)
s.x <- function(x,theta) theta[4]/(2*sqrt(theta[3]+theta[4]*x))
s.xx <- function(x,theta) -theta[4]^2/(4*(theta[3]+theta[4]*x)^1.5)
# AIC for competing model 1
sdeAIC(X, NULL, b, s, b.x, s.x, s.xx, guess=c(1,1,1,1),
lower=rep(1e-3,4), method="L-BFGS-B")
s <- function(x,theta) (theta[3]+theta[4]*x)^theta[5]
s.x <- function(x,theta)
theta[4]*theta[5]*(theta[3]+theta[4]*x)^(-1+theta[5])
s.xx <- function(x,theta) (theta[4]^2*theta[5]*(theta[5]-1)
*(theta[3]+theta[4]*x)^(-2+theta[5]))
# AIC for competing model 2
sdeAIC(X, NULL, b, s, b.x, s.x, s.xx, guess=c(1,1,1,1,1),
lower=rep(1e-3,5), method="L-BFGS-B")

set.seed(123)
theta0 <- c(0.89218*0.09045,0.89218,sqrt(0.032742))
theta1 <- c(0.89218*0.09045/2,0.89218,sqrt(0.032742/2))
# we test the true model against two competing models
b <- function(x,theta) theta[1]-theta[2]*x
b.x <- function(x,theta) -theta[2]
s <- function(x,theta) theta[3]*sqrt(x)
s.x <- function(x,theta) theta[3]/(2*sqrt(x))
s.xx <- function(x,theta) -theta[3]/(4*x^1.5)
X <- sde.sim(X0=rsCIR(1, theta1), N=1000, delta=1e-3, model="CIR", theta=theta1)
sdeDiv(X=X, theta0 = theta0, b=b, s=s, b.x=b.x, s.x=s.x, s.xx=s.xx, method="L-BFGS-B",
lower=rep(1e-3,3), guess=c(1,1,1))
sdeDiv(X=X, theta0 = theta1, b=b, s=s, b.x=b.x, s.x=s.x, s.xx=s.xx, method="L-BFGS-B",
lower=rep(1e-3,3), guess=c(1,1,1))
lambda <- -1.75
myphi <- expression( (x^(lambda+1) -x - lambda*(x-1))/(lambda*(lambda+1)) )
sdeDiv(X=X, theta0 = theta0, phi = myphi, b=b, s=s, b.x=b.x, s.x=s.x, s.xx=s.xx, method="L-BFGS-B",
lower=rep(1e-3,3), guess=c(1,1,1))
sdeDiv(X=X, theta0 = theta1, phi = myphi, b=b, s=s, b.x=b.x, s.x=s.x, s.xx=s.xx, method="L-BFGS-B",
lower=rep(1e-3,3), guess=c(1,1,1))

## Not run:
set.seed(123)
d <- expression(-1*x)
s <- expression(2)
sde.sim(drift=d, sigma=s,N=50,delta=0.01) -> X
S <- function(t, x, theta) sqrt(theta[2])
B <- function(t, x, theta) -theta[1]*x
true.loglik <- function(theta) {
DELTA <- deltat(X)
lik <- 0
for(i in 2:length(X))
lik <- lik + dnorm(X[i], mean=X[i-1]*exp(-theta[1]*DELTA),
sd = sqrt((1-exp(-2*theta[1]*DELTA))*
theta[2]/(2*theta[1])),TRUE)
lik
}
xx <- seq(-10,10,length=20)
sapply(xx, function(x) true.loglik(c(x,4))) -> py
sapply(xx, function(x) EULERloglik(X,c(x,4),B,S)) -> pz
sapply(xx, function(x) SIMloglik(X,c(x,4),B,S,M=10000,N=5)) -> pw
plot(xx,py,type="l",xlab=expression(beta),
ylab="log-likelihood",ylim=c(0,15)) # true
lines(xx,pz, lty=2) # Euler
lines(xx,pw, lty=3) # Simulated
## End(Not run)

set.seed(123);
# Kessler’s estimator for O-H process
K.est <- function(x) {
n.obs <- length(x)
n.obs/(2*(sum(x^2)))
}
# Least squares estimators for the O-H process
LS.est <- function(x) {
n <- length(x) -1
k.sum <- sum(x[1:n]*x[2:(n+1)])
dt <- deltat(x)
ifelse(k.sum>0, -log(k.sum/sum(x[1:n]^2))/dt, NA)
}
d <- expression(-1 * x)
s <- expression(1)
x0 <- rnorm(1,sd=sqrt(1/2))
sde.sim(X0=x0,drift=d, sigma=s,N=2500,delta=0.1) -> X
# Kessler’s estimator as estimating function
f <- list(expression(2*theta*x^2-1))
simple.ef(X, f, lower=0, upper=Inf)
K.est(X)
# Least Squares estimator as estimating function
f <- list(expression(x*(y-x*exp(-0.1*theta))))
simple.ef(X, f, lower=0, upper=Inf)
LS.est(X)

set.seed(123)
d <- expression(10 - x)
s <- expression(sqrt(x))
x0 <- 10
sde.sim(X0=x0,drift=d, sigma=s,N=1500,delta=0.1) -> X
# rather difficult problem unless a good initial guess is given
d <- expression(alpha + theta*x)
s <- expression(x^gamma)
h <- list(expression(x), expression(x^2), expression(x^2))
simple.ef2(X, d, s, h, lower=c(0,-Inf,0), upper=c(Inf,0,1))


# ぱらぱらめくる『Simulation and Inference for Stochastic Differential Equations: With R Examples』1. 確率過程とそれを理解するための基礎。Ito過程。確率微分方程式とは

• 1. Stochastic Processes and Stochastic Differential Equations
• 確率空間と確率過程の数学的記法
• 確率変数の平均・分散・モーメント
• シミュレーションのための道具
• 疑似乱数列
• モンテカルロ法
• 分散を小さくするための手法
• モンテカルロはよい方法だが、あまりにナイーブにやりすぎると、収束が遅くて使えない(使いにくい)。その理由は、分散(信頼区間)が広すぎることに起因する(ことが多い)。それを解決するために、「早く収束させる」技術がいくつかある
• 確率過程は「経路・道」
• 確率過程のばらつき
• 確率過程の特徴づけ
• モーメント、共分散、増分
• ブラウン運動
• ３つのシミュレーション方法
• 増分を正規分布でとる
• 酔歩(-1,1)の繰り返し
• Karhunen-Loeve 拡張(たくさんのものの足し合わせの極限としての定義がある。その足し合わせをある程度すれば、シミュレーションとしては精度が出る、といった感じの「定義とシミュレーション法」)
• どれがよい？
• ある時点の状態分布を取るだけなら、増分を正規分布でとるだけで十分
• 経過を問題にするならKarhunen-Loeve 拡張がよい
• ３列はそれぞれ、左から「増分＝正規分布」「酔歩」「Karhunen-Loeve拡張」
• ３行はそれぞれ、上から、時刻数が100,1000,100、Karhunen-Loeve拡張での重ね合わせ数50,50,10000
• 時刻数を増やすと経過時間が伸びるのが「増分・正規」「酔歩」。「Karhunen-Loeve拡張」は増えていない
• 「Karhunen-Loeve拡張」は時刻数を増やしても、重ね合わせ数を増やしても「細かさが増える」

N1<-100
N2<-1000
T<-1
Delta1<-T/N1
Delta2<-T/N2
W1.1<-W2.1<-W3.1<-numeric(N1+1)
W1.2<-W2.2<-W3.2<-numeric(N2+1)
W3.3<-numeric(N1+1)
t1<-seq(from=0,to=T,length=N1+1)
t2<-seq(from=0,to=T,length=N2+1)
# 単位時間当たりの増分を正規分布から乱数発生する
set.seed(123)
W1.1<-c(0,cumsum(rnorm(N1)*sqrt(Delta)))
set.seed(123)
W1.2<-c(0,cumsum(rnorm(N2)*sqrt(Delta)))

# 酔歩
set.seed(123)
W2.1<-c(0,cumsum(sample(c(-1,1),N1,replace=TRUE)))/sqrt(N1)
set.seed(123)
W2.2<-c(0,cumsum(sample(c(-1,1),N2,replace=TRUE)))/sqrt(N2)

# Karhunen-Loeve expansion
# 足し合わせる分布数を指定するパラメタが増える
n1<-50
n2<-10000
set.seed(123)
Z1<-rnorm(n1)
phi<-function(i,t,T){
(2*sqrt(2*T))/((2*i+1)*pi) * sin(((2*i+1)*pi*t)/(2*T))
}
for(i in 1:n1){
W3.1<-W3.1+Z1[i]*phi(i,t1,T)
}
set.seed(123)
Z2<-rnorm(n2)
for(i in 1:n2){
W3.2<-W3.2+Z2[i]*phi(i,t2,T)
}
for(i in 1:n2){
W3.3<-W3.3+Z2[i]*phi(i,t1,T)
}

ylim<-range(W1.1,W2.1,W3.1,W1.2,W2.2,W3.2,W3.3)
par(mfcol=c(3,3))
plot(t1,W1.1,type="l",ylim=ylim)
plot(t2,W1.2,type="l",ylim=ylim)
plot(t1,W1.1,type="l",ylim=ylim)
plot(t1,W2.1,type="l",ylim=ylim)
plot(t2,W2.2,type="l",ylim=ylim)
plot(t1,W2.1,type="l",ylim=ylim)
plot(t1,W3.1,type="l",ylim=ylim)
plot(t2,W3.2,type="l",ylim=ylim)
plot(t1,W3.3,type="l",ylim=ylim)
par(mfcol=c(1,1))

• いたるところ微分不可能
• ブラウン運動の拡張例
• 確率的積分と確率的微分方程式
• いたるところ微分不可能
• 微分の式はできるが(確率的微分方程式)、そこに、無限小時間で発散する項が入っているので、意味がない…。意味がないけれど、微分方程式の係数は求めたい…
• 確率微分方程式の例：がWiener過程になっていて、これは、「いたるところ微分不可能」なので、この式自体をそのまま扱うことはしない
• 積分はできる過程はある(Ito sum, Ito integralが可能な過程: Ito 過程)
• 部分ごとに近似した積分を取る
• 部分をどんどん狭くして、極限を取る
• Ito 過程の特徴
• 拡散過程
• ドリフト(どこかへ向かっていく項)と拡散の２項からなる
• いくつかの過程(取扱いを可能にするための)
• 不等式で示されるいくつかの特徴
• エルゴード的(Wiki)
• マルコフ性(Wiki)
• 拡散過程のinfinitesimal generator
• Ito 公式
• 確率過程の計算に重要
• シミュレーションにも重要
• テイラー展開の確率過程版
• Ito 公式はよいが、実用(シミュレーションにも、微分方程式の解法にも(?))には工夫が必要
• 変換する
• 測度を変える
• 尤度を計算しやすくする
• Parametric families of stochastic processes(使いやすいいくつかの確率過程)
• The Ornstein-Uhlenbech or Vasicek process
• The Black-Scholes-Merton or geometric Brownian motion model
• The Cox-Ingersoll-Ross model
• The Chan-Karolyi-Longstaff-Sanders (CKLS) family of models
• The modified CIR and hyperbolic processes
• The hyperbolic processe
• The nonlinear mean reversion Ait-Sahaliamodel
• Double-well potential
• The Jacobi diffusion process
• Ahn and Gao model or inverse of Feller's square root model
• Rarial Ornstein-Uhlenbeck process
• Pearson diffusions
• Another classification of linear stochastic systems
• One epidemic model
• The stochastic cusp catastrophe model
• Exponential families of diffusions
• Generalized inverse gaussian diffusions

# ぱらぱらめくる『Simulation and Inference for Stochastic Differential Equations: With R Examples』2. 確率微分方程式のシミュレーション

• 2. Numerical Nethods for SDE
• 離散的増分でシミュレーションする
• 係数から、連続変化の(任意時刻〜離散時刻)の値をシミュレーションする
• sde.sim()関数はシミュレーション関数
• method引数の"euler","milstein","KPS","milstein2"は離散的増分シミュレーション、"cdist","ozaki","shoji","EA"は係数からの連続変化シミュレーション
• を引数 drift, sigmaとして指定するのが基本
• drift, sigmaには式が入れられる

# ぱらぱらめくる『Simulation and Inference for Stochastic Differential Equations: With R Examples』3. 確率過程への確率微分方程式モデルの当て嵌めとそのパラメタ推定。適切なモデルの選択

• 3. Parametric Estimation
• 推定のためには、推定結果を評価する方法が必要。それが「尤度」関係
• Likelihoodの計算
• Exact likelihoodを計算
• Pseudo-likelihoodを計算
• Likelihoodの近似計算
• 観察データの何を基に推定するか
• 変化の動き自体？、時間と分散の関係？…
• 推定の仕方は色々な発表があり、それに対して、関数を作ってくれている
• simple.ef()
• simple.ef2()
• linear.mart.ef()
• gmm()