Numerical Solution of Stochastic Differential Equations, Peter E.Kloeden Eckhard Platen, Springer-Verlag

```{r}

#最尤推定から生成作用素を計算する　p.33

n=3
t=10
num=array(0,dim=c(n^2,2))
count=1

for(i in 1:n){
for(j in 1:n){  
num[count,]=c(i,j) 
count=count+1  
}}
Y=matrix(sample(c(1:9),t*n^2,replace=T),ncol=t)
Y=t(t(Y)/apply(Y,2,sum))


#Pij(t)=(1-exp(-t*lam^2))/(1+sum(1-exp(-t*lam^2)))
#Pii(t)=1/(1+sum(1-exp(-t*lam^2)))
#loglik=sum(Yijt*log(Pij(t)))

loglik=function(lam){
logliklihood=array(0,dim=c(n,n))  
for(k in 1:length(t)){  
P=(1-exp(-t[k]*lam^2))/(1+apply(1-exp(-t[k]*lam^2),2,sum))  
diag(P)=1/(1+apply(1-exp(-t[k]*lam^2),2,sum)) 
logliklihood=logliklihood+matrix(Y[,1],ncol=n)*log(P)
}
return(sum(logliklihood))
}  
  
lam=array(0.1,dim=c(n,n))
diag(lam)=0

ite=10^5
eta=0.01
h=0.001

for(l in 1:ite){
coef=lam  
for(i in 1:n){
for(j in 1:n){  
if(i!=j){
coef_sub=coef;coef_sub[i,j]=coef_sub[i,j]+h
lam[i,j]=lam[i,j]+eta*(loglik(coef_sub)-loglik(coef))/h
}}}
print(loglik(lam))
}

#生成作用素を計算する
A=lam^2
diag(A)=-apply(A,1,sum)



```

```{r}


#stratonovich integral p.157

h=function(s,x){
return(exp(-(s*x)^2))  
}

#stochastic differential equation
#dxt=a*dt+b*dwt

a=function(s,x){
return(-2*x+s)  
}

b=function(s,x){
return(-3*s+x)  
}


dt=0.001
t=seq(dt,1,dt)

ite=1000
mesh=0.0001

Ito=c()
stratonovich=c()

for(k in 1:ite){

X=0
W=rnorm(1,0,sqrt(dt))
stratonovich_integral=0
Ito_integral=0

for(l in 1:length(t)){
W_pre=W
W=rnorm(1,0,sqrt(dt))+W_pre
X_pre=X
abar=a(t[l]-dt,X_pre)-b(t[l]-dt,X_pre)*(b(t[l]-dt,X_pre+mesh)-b(t[l]-dt,X_pre))/(2*mesh)
X=X_pre+abar*(dt)+b(t[l]-dt,X_pre)*(W-W_pre)+b(t[l]-dt,X_pre)*(b(t[l]-dt,X_pre+mesh)-b(t[l]-dt,X_pre))/(2*mesh)
stratonovich_integral=stratonovich_integral+h(t[l]-dt,(X+X_pre)/2)*(W-W_pre)
X=X_pre+a(t[l]-dt,X_pre)*(dt)+b(t[l]-dt,X_pre)*(W-W_pre)
Ito_integral=Ito_integral+h(t[l]-dt,X_pre)*(W-W_pre)
}

Ito=c(Ito,Ito_integral)

stratonovich=c(stratonovich,(stratonovich_integral))

#print(k)
#print(mean(Ito))
#print(mean(stratonovich))

}

plot(stratonovich,main=paste0("ave:",mean(stratonovich)))
plot(Ito,main=paste0("ave:",mean(Ito)))



```

```{r}

#Turbulent Diffusion and Radio-Astronomy p.259
#one dimension

dt=0.001
T=3
t=seq(dt,T,dt)
X=0
V=0
sig=10
W=0
X_vec=c()
V_vec=c()

for(l in 1:length(t)){

W_pre=W
W=rnorm(1,0,dt)+W_pre
X=X+V*dt
V=V-V*dt/T+sig*(W-W_pre)

X_vec=c(X_vec,X)
V_vec=c(V_vec,V)
}

plot(t,X_vec,xlim=c(min(t),max(t)),ylim=c(min(c(V_vec,X_vec)),max(c(V_vec,X_vec))),col=2,xlab="t",ylab="values")

par(new=T)

plot(t,V_vec,xlim=c(min(t),max(t)),ylim=c(min(c(V_vec,X_vec)),max(c(V_vec,X_vec))),col=3,xlab="t",ylab="values")



```

```{r}

#p.308

#Euler approximation

#dxt=a*Xt*dt+b*Xt*dWt
#厳密解Xt=X0*exp((a-b^2/2)*t+b*Wt)

X=0.1
Y0=X
W=0
dt=0.01
t=seq(dt,10,dt)
a=1
b=2

X_path=c()
Y_path=c()

for(l in 1:length(t)){
W_pre=W  
W=W+rnorm(1,0,sqrt(dt))  
X=X+a*X*dt+b*X*(W-W_pre)  
Y=Y0*exp((a-b^2/2)*t[l]+b*W)
X_path=c(X_path,X)
Y_path=c(Y_path,Y)
}

#差分による近似解
plot(t,X_path,xlim=c(min(t),max(t)),ylim=c(min(c(Y_path,X_path)),max(c(Y_path,X_path))),col=2,xlab="t",ylab="value")

par(new=T)

#厳密解
plot(t,Y_path,xlim=c(min(t),max(t)),ylim=c(min(c(Y_path,X_path)),max(c(Y_path,X_path))),col=3,xlab="t",ylab="value")

```