金融モデルにおける推定と最適化　乾孝治・室町幸雄先生　朝倉書店

```r

#p.44 GMM(Generalized Moment Method:一般化モーメント法)
#(2.43)を書き換えたYi=(sita1-sita2^2/2)+sita2*epでのパラメータ推定
Y=(20-4^2/2)+rnorm(300,0,4)

cost=function(sita){
h=rbind(Y-(sita[1]-sita[2]^2/2),Y^2-(sita[2]+(sita[1]-sita[2]^2/2)^2),Y^3-3*(sita[1]-sita[2]^2/2)*sita[2]-(sita[1]-sita[2]^2/2)^3)
S=diag(1,3)
for(j in 1:ncol(h)){
S=S+t(t(c(h[,j])))%*%t(c(h[,j]))  
}
S=S-diag(1,3)
S=S/ncol(h)
W=solve(S)
g=apply(h,1,mean)
Jt=as.numeric(g%*%W%*%g)
return(Jt)
}



Sita=c(mean(Y),sd(Y))
h=0.001
learning_rate=10^(-1)
ite=5000

for(l in 1:ite){
dSita_vec=c()
for(j in 1:length(Sita)){
Sita_sub=Sita
Sita_sub[j]=Sita_sub[j]+h
dSita_vec=c(dSita_vec,(cost(Sita_sub)-cost(Sita))/h)
}
m=-learning_rate*dSita_vec
Sita=Sita+m
print(cost(Sita))
}

predict=Sita[1]-Sita[2]^2/2+rnorm(length(Y),0,Sita[2]^2)

plot(c(1:length(Y)),predict,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),xlab="index",ylab="objective variables",main="stochastic differential equation",col=2)

par(new=T)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),xlab="index",ylab="objective variables",main="stochastic differential equation",col=3)


```

```r

#逐次二次計画法 p.89

f=function(x){
return(2*x[1]^2+3*x[2]^2+5)  
}

A=matrix(c(4,3,1,1,2,2,-1,0,0,0,-1,0,0,0,-1),nrow=3)
b=c(13,16,8)
cx=function(x){
return(A%*%x-b)
}

L=function(x,v){
return(f(x)-sum(cx(x)*v))
}


ite=100
X=rep(1,5)
V=rep(1,3)
h=10^(-4)
B=diag(1,length(X))
eta=10^(-2)

BFGS=0
DFP=0
SROU=1
SM=0

for(l in 1:ite){
Xpre=X
  
y=c();dcx=array(0,dim=c(length(X),length(V)));df=c()
for(j in 1:length(X)){
Xsub=X;Xsub[j]=Xsub[j]+h
y=c(y,-(L(Xsub,V)-L(X,V))/h)
dcx[j,]=((cx(Xsub)-cx(X))/h)
df=c(df,(f(Xsub)-f(X))/h)
}

H=(B)
matrix=rbind(cbind(H,-dcx),cbind(t(dcx),array(0,dim=c(length(V),length(V)))))
vector=-c(df,cx(X))
solve=solve(matrix)%*%vector
d=solve[1:length(X)]
V=V+eta*solve[-c(1:length(X))]
X=X+eta*d
s=X-Xpre
  
for(j in 1:length(X)){
Xsub=X;Xsub[j]=Xsub[j]+h
y[j]=y[j]+(L(Xsub,V)-L(X,V))/h
}

if(BFGS==1){

#p.90 参照
if(sum(s*y)>=0.2*as.numeric(c(s)%*%B%*%s)){
sita=1  
}else{
sita=0.8*as.numeric(c(s)%*%B%*%s)/(as.numeric(c(s)%*%B%*%s)-sum(s*y))  
}

y=sita*y+(1-sita)*B%*%s
B=B-as.numeric(B%*%s%*%t(B%*%s))/as.numeric(c(s)%*%B%*%s)+y%*%t(y)/sum(s*y)

}

if(DFP==1){
B=B-((B%*%s)%*%t(y)+y%*%t(B%*%s))/sum(s*y)+(1+as.numeric(c(s)%*%B%*%c(s))/sum(s*y))*y%*%t(y)/sum(s*y)
}

if(SROU==1){
B=B+(y-B%*%s)%*%t(y-B%*%s)/sum(s*c(y-B%*%s))  
}

if(SM==1){
H=solve(B)
H=H-(H%*%y%*%t(s)+s%*%t(H%*%y))/sum(s*y)+(1+as.numeric(c(y)%*%H%*%c(y))/sum(s*y))*s%*%t(s)/sum(s*y)
B=solve(H)
}

print(L(X,V))
#print(sum(cx(X)))
}




```

```r


#p.116 最尤法
#オイラー法の離散化より確率微分方程式のパラメータ推定
data("AirPassengers")
r=AirPassengers
r_lag=r[1:(length(r)-1)]
r_pred=r[2:length(r)]

cost=function(alpha,beta,gam,sigma){
pred=alpha+(beta+1)*r_lag 
ep=r_pred-pred  
loglike=-length(r_lag)*log(2*pi*sigma^2)/2-gam*sum(log(r_lag))-sum(ep^2/(2*(sigma*r_lag^gam)^2))
return(-loglike)
}

param=solve(matrix(c(sum(r_lag),length(r_pred),-sum(r_lag^2),-sum(r_lag)),nrow=2))%*%c(sum(r_lag^2)-sum(r_lag*r_pred),sum(r_lag)-sum(r_pred))

a=param[1]
b=param[2]
gamma=1.1
sig=1
h=0.001
learning_rate=10^(-6)
ite=10^4

for(l in 1:ite){
 
da=(cost(a+h,b,gamma,sig)-cost(a,b,gamma,sig))/h
db=(cost(a,b+h,gamma,sig)-cost(a,b,gamma,sig))/h
dgamma=(cost(a,b,gamma+h,sig)-cost(a,b,gamma,sig))/h
dsig=(cost(a,b,gamma,sig+h)-cost(a,b,gamma,sig))/h

a=a-learning_rate*da
b=b-learning_rate*db
gamma=gamma-learning_rate*dgamma
sig=sig-learning_rate*dsig

print(cost(a,b,gamma,sig))
}

#E[rt] ep_tが平均0の場合
predict=a+(b+1)*r_lag

plot(c(1:length(r_lag)),predict,xlim=c(1,length(r_lag)),ylim=c(min(r_pred),max(r_pred)),xlab="index",ylab="objective variables",main="stochastic differential equation",col=2)

par(new=T)

plot(c(1:length(r_lag)),r_pred,xlim=c(1,length(r_lag)),ylim=c(min(r_pred),max(r_pred)),xlab="index",ylab="objective variables",main="stochastic differential equation",col=3)



```

```r

#p.116 GMM(Generalized Moment Method:一般化モーメント法)
#オイラー法の離散化より確率微分方程式のパラメータ推定
data("AirPassengers")
r=AirPassengers
r_lag=r[1:(length(r)-1)]
r_pred=r[2:length(r)]

cost=function(alpha,beta,gam,sigma){
pred=alpha+(beta+1)*r_lag 
ep=r_pred-pred  
eta=ep^2-(sigma*r_lag^gam)^2  
h=rbind(ep,ep*r_lag,eta,eta*r_lag)
S=diag(1,4)
for(j in 1:ncol(h)){
S=S+t(t(c(h[,j])))%*%t(c(h[,j]))  
}
S=S-diag(1,4)
S=S/ncol(h)
W=solve(S)
g=apply(h,1,mean)
Jt=as.numeric(g%*%W%*%g)
return(Jt)
}

a=1
b=1
gamma=0.2
sig=1
h=0.001
learning_rate=10^(-2)
ite=10^4

for(l in 1:ite){
 
da=(cost(a+h,b,gamma,sig)-cost(a,b,gamma,sig))/h
db=(cost(a,b+h,gamma,sig)-cost(a,b,gamma,sig))/h
#dgamma=(cost(a,b,gamma+h,sig)-cost(a,b,gamma,sig))/h
dsig=(cost(a,b,gamma,sig+h)-cost(a,b,gamma,sig))/h

a=a-learning_rate*da
b=b-learning_rate*db
#gamma=gamma-learning_rate*dgamma
sig=sig-learning_rate*dsig

#print(cost(a,b,gamma,sig))
}

#E[rt] ep_tが平均0の場合
predict=a+(b+1)*r_lag

plot(c(1:length(r_lag)),predict,xlim=c(1,length(r_lag)),ylim=c(min(r_pred),max(r_pred)),xlab="index",ylab="objective variables",main="stochastic differential equation",col=2)

par(new=T)

plot(c(1:length(r_lag)),r_pred,xlim=c(1,length(r_lag)),ylim=c(min(r_pred),max(r_pred)),xlab="index",ylab="objective variables",main="stochastic differential equation",col=3)


```

```r

#p.132 3次スプラインにおける重回帰分析

data("AirPassengers")
dat=data.frame(AirPassengers)
N=130
Pm=dat[1:N,]
c=1
M=100
t=c(0:M)

fjt=function(time){
f=c()
for(j in 2:(length(t)-1)){
if(t[j-1]<time & time<=t[j]){
f=c(f,(time-t[j-1])^3/(6*(t[j]-t[j-1])))  
}
if(t[j]<time & time<=t[j+1]){
c=t[j]-t[j-1]
e=time-t[j]
f=c(f,c^2/6+c*e/2+e^2/2-e^3/(6*(t[j+1]-t[j])))
}
if(t[j+1]<time){
f=c(f,(t[j+1]-t[j-1])*((2*t[j+1]-t[j]-t[j-1])/6+(time-t[j+1])/2))  
}
if(time<=t[j-1]){f=c(f,0)}
}
return(f)
}


tvec=seq(1,N,1)
N=length(tvec)

fmat=array(0,dim=c(N,M-1))
for(j in 1:N){
fmat[j,]=fjt(tvec[j])  
}


a=solve(t(fmat)%*%fmat)%*%t(fmat)%*%Pm
cor(Pm,fmat%*%a)

predict=fmat%*%a
Y=Pm

plot(c(1:length(Y)),predict,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),xlab="index",ylab="objective variables",main="Regression",col=2)
par(new=T)
plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(Y),max(Y)),xlab="index",ylab="objective variables",main="Regression",col=3)


```

```r

#p.134 Vasicek and Fong
#作成中

data("AirPassengers")
dat=data.frame(AirPassengers)
N=50
Y=dat[1:N,]
Y=cumsum(Y);Y=Y/max(Y)
c=1/N
t=seq(c,1,c)

fjt=function(time){
f=c()
for(j in 2:(length(t)-1)){
if(t[j-1]<time & time<=t[j]){
f=c(f,(time-t[j-1])^3/(6*(t[j]-t[j-1])))  
}
if(t[j]<time & time<=t[j+1]){
c=t[j]-t[j-1]
e=time-t[j]
f=c(f,c^2/6+c*e/2+e^2/2-e^3/(6*(t[j+1]-t[j])))
}
if(t[j+1]<time){
f=c(f,(t[j+1]-t[j-1])*((2*t[j+1]-t[j]-t[j-1])/6+(time-t[j+1])/2))  
}
if(time<=t[j-1]){f=c(f,0)}
}
return(f)
}


f=function(W,s){
alpha=W[1];beta=W[-1]
tvec=seq(c,1,c)
fmat=array(0,dim=c(N,N-2))
for(j in 1:N){
fmat[j,]=fjt(1-exp(-alpha*tvec[j]))  
}
if(s==1){
fs=fmat%*%beta
}else{
fs=solve(t(fmat)%*%fmat+diag(1,ncol(fmat)))%*%t(fmat)%*%Y  
}
return(fs)  
}

sita=c(1,f(sita,2))
ite=10^3
eta=10^(-1)
lam=1
h=0.001

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub,1)-f(vec,1))/h
Z[,j]=Z_sub
}
b=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita,1))
sita=c(sita+eta*b)
print(sum((Y-f(sita,1))^2))
}

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(c(Y,f(sita,1))),max(c(Y,f(sita,1)))),xlab="Index",ylab="Objective Variable",col=2,type="o")
par(new=T)
plot(c(1:length(Y)),f(sita,1),xlim=c(1,length(Y)),ylim=c(min(c(Y,f(sita,1))),max(c(Y,f(sita,1)))),xlab="Index",ylab="Objective Variable",col=3,type="o")



```