ぱらぱらめくる『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))