保険と金融の数理　共立出版　室井芳史先生

```r

#p.54

#オイラー法による遷移確率の計算
#0:死亡,1:入院,2:健康

lam1=0
lam2=0
lam10=0
lam20=0
lam12=0
lam21=0


b1=25.5
b2=15.5
b10=150
b20=225
b12=20.5
b21=12.5

m1=86.34
m2=86.34
m10=86.34
m20=86.34
m12=86.34
m21=86.34

mu12=function(x){return(lam12+exp((x-m12)/b12)/b12)}
mu21=function(x){return(lam21+exp((x-m21)/b21)/b21)}
mu10=function(x){return(lam10+exp((x-m10)/b10)/b10)}
mu20=function(x){return(lam20+exp((x-m20)/b20)/b20)}
mu1=function(x){return(mu12(x)+mu10(x))}
mu2=function(x){return(mu21(x)+mu20(x))}


#初期確率
p10=0
p20=0
p11=1
p12=0
p21=0
p22=1

ite=1200
dt=0.1
t=0

t_vec=c(t)
p22_vec=c(p22)
p21_vec=c(p21)
p20_vec=c(p20)
p12_vec=c(p12)
p11_vec=c(p11)
p10_vec=c(p10)

for(l in 1:ite){
  
t=t+dt;t_vec=c(t_vec,t)  
p11=p11+dt*(p12*mu21(t)-p11*mu1(t));p11_vec=c(p11_vec,ifelse(p11>0,p11,ifelse(p11<1,p11,1)))
p22=p22+dt*(-p22*mu2(t)+p21*mu12(t));p22_vec=c(p22_vec,ifelse(p22>0,p22,ifelse(p22<1,p22,1)))
p21=p21+dt*(p22*mu21(t)-p21*mu1(t));p21_vec=c(p21_vec,ifelse(p21>0,p21,ifelse(p21<1,p21,1)))
p20=p20+dt*(p22*mu20(t)+p21*mu10(t));p20_vec=c(p20_vec,ifelse(p20>0,p20,ifelse(p20<1,p20,1)))
p12=p12+dt*(-p12*mu2(t)+p11*mu12(t));p12_vec=c(p12_vec,ifelse(p12>0,p12,ifelse(p12<1,p12,1)))
p10=p10+dt*(p12*mu20(t)+p11*mu10(t));p10_vec=c(p10_vec,ifelse(p10>0,p10,ifelse(p10<1,p10,1)))
}

plot(t_vec,p11_vec,xlab="t",ylab="p",xlim=c(min(t_vec),max(t_vec)),main="p11_vec",col=2)

plot(t_vec,p12_vec,xlab="t",ylab="p",xlim=c(min(t_vec),max(t_vec)),main="p12_vec",col=2)

plot(t_vec,p22_vec,xlab="t",ylab="p",xlim=c(min(t_vec),max(t_vec)),main="p22_vec",col=2)

plot(t_vec,p21_vec,xlab="t",ylab="p",xlim=c(min(t_vec),max(t_vec)),main="p21_vec",col=2)

plot(t_vec,p10_vec,xlab="t",ylab="p",xlim=c(min(t_vec),max(t_vec)),main="p10_vec",col=2)

plot(t_vec,p20_vec,xlab="t",ylab="p",xlim=c(min(t_vec),max(t_vec)),main="p20_vec",col=2)


```

```r

#p.59

m=86.34
b=9.5
lam=0

mux=function(x){lam+exp((x-m)/b)/b}

t=0
p=1
dt=0.1
ite=1000

t_vec=c(t)
p_vec=c(p)

for(l in 1:ite){

t=t+dt  
t_vec=c(t_vec,t)

s1=-mux(t)
s2=-mux(t+dt/2)*(p+dt*s1/2)
s3=-mux(t+dt/2)*(p+s2*dt/2)
s4=-mux(t+dt)*(p+s3*dt)
p=p+dt*(s1+2*s2+2*s3+s4)/6
p_vec=c(p_vec,p)
  
}

plot(t_vec,p_vec,col=2,ylab="p")


```

```r


#確率微分方程式のパラメータ推定 p.106
#オイラー法

data("AirPassengers")

x=log(c(AirPassengers))
x0=100
y=c(x0,x[1:(length(x)-1)])

p=2

mu=function(a,X){
vec=c(1,X^c(1:p))
return((sum(a*vec)))
}

A=rep(0.0001,p+1)
sigma=1

ite=10^4
eta=10^(-6)
dt=0.1
t=0
h=0.001


for(l in 1:ite){
  
muvec=c()
dmuvec=array(0,dim=c(length(x),length(A)))

for(i in 1:length(y)){
muvec=c(muvec,mu(A,y[i]))  
for(j in 1:length(A)){
Asub=A;Asub[j]=Asub[j]+h
dmuvec[i,j]=(mu(Asub,y[i])-mu(A,y[i]))/h
}}

dsig=-pi*dt*2*sigma*length(y)+sum(c(c(x-y-muvec*dt)^2/c(dt*sigma^3)))  
dmu=t(dmuvec)%*%c(dt*c(x-y-muvec*dt)/c(dt*sigma^2))

A=A+eta*(dmu)
sigma=sigma+eta*(dsig)

loglik=-pi*dt*sum(sigma^2)-sum((x-y-muvec*dt)^2/(2*dt*sigma^2))
print(loglik)
}



plot(c(1:length(y)),y+muvec*dt+rnorm(length(y),sd=sqrt(dt*sigma^2)),xlim=c(1,length(y)),ylim=c(min(c(x,y+muvec*dt))-3,max(c(x,y+muvec*dt))+3),xlab="samples",ylab="values",col=2)

par(new=T)

plot(c(1:length(y)),x,xlim=c(1,length(y)),ylim=c(min(c(x,y+muvec*dt))-3,max(c(x,y+muvec*dt))+3),xlab="samples",ylab="values",col=3)


```

```r
#確率微分方程式のパラメータ推定
#局所線形化(3.56)式　p.107
#x=a0+a1*y+a2*y^2+a3*t

data("AirPassengers")

x=log(c(AirPassengers))
x0=log(100)
y=c(x0,x[1:(length(x)-1)])

h=0.001
p=2
dt=0.1
n=length(y)
A=rep(0.01,p+2)
sig=1

mu=function(a,yy,tt){
Y=array(0,dim=c(length(yy),p+1))
for(j in 1:(p+1)){
Y[,j]=yy^(j-1)  
}
return(Y%*%a[1:(p+1)]+a[-c(1:(p+1))]*tt)  
}

loglik=function(a,sigma,yy){
Mn=(sigma^2)*(mu(a,yy+h,dt*(c(1:n)-1))-2*mu(a,yy,dt*(c(1:n)-1))+mu(a,yy-h,dt*(c(1:n)-1)))/(2*h^2)+(mu(a,yy,dt*(c(1:n)-1)+h)-mu(a,yy,dt*(c(1:n)-1)))/h 
Ln=(mu(a,yy+h,dt*(c(1:n)-1))-mu(a,yy,dt*(c(1:n)-1)))/h  
En=yy+mu(a,yy,dt*(c(1:n)-1))*(exp(Ln*dt)-1)/Ln+Mn*((exp(Ln*dt)-1)-Ln*dt)/(Ln^2) 
Vn=(exp(2*Ln*dt)-1)*sigma^2/(2*Ln) 
return(-sum((x-En)^2/Vn)/2-sum(log(2*pi*(Vn+1))/2))
}


ite=10^4
eta=0.0001
for(l in 1:ite){
for(j in 1:length(A)){
Asub=A;Asub[j]=Asub[j]+h
A[j]=A[j]+eta*(loglik(Asub,sig,y)-loglik(A,sig,y))/h
}  
sig=sig+eta*(loglik(A,sig+h,y)-loglik(A,sig,y))/h    
print(loglik(A,sig,y))
}

muvec=mu(A,y,dt*(c(1:n)-1))

plot(c(1:length(y)),y+muvec*dt+rnorm(length(y),sd=sqrt(dt*sig^2)),xlim=c(1,length(y)),ylim=c(min(c(x,y+muvec*dt))-3,max(c(x,y+muvec*dt))+3),xlab="samples",ylab="values",col=2)

par(new=T)

plot(c(1:length(y)),x,xlim=c(1,length(y)),ylim=c(min(c(x,y+muvec*dt))-3,max(c(x,y+muvec*dt))+3),xlab="samples",ylab="values",col=3)



```