多変量ノンパラメトリック回帰と視覚化

#ARCH　p.142

data("AirPassengers")

t=length(AirPassengers)

p=4

Y=AirPassengers

f=function(param){

cost=0  

val=rep(0,length(param))

for(j in 1:(t-p+1)){  

y=c(1,Y[j:(j+p-1)]^2)

sig=sum(y*param)

cost=cost+log(sig)+Y[j+p-1]^2/(sig)

#for(k in 1:length(param)){

#if(k==1){

#val[k]=val[k]+1/param[k]+Y[j+p-1]^2/sig  

#}else{

#val[k]=val[k]+Y[j+k-1]^2+(Y[j+p-1]*Y[j+k-1])^2/sig  

#}

#}

}

return(cost)

}

ite=10^3

eta=10^(-5)

h=0.001

params=rep(1,p+1)

df_pre=rep(1,length(params))

for(l in 1:ite){

df=c()  

for(j in 1:length(params)){

vec=params;vec[j]=vec[j]+h  

df=c(df,(f(vec)-f(params))/h)

params[j]=params[j]-eta*(f(vec)-f(params))/h

}  

eta=1/(sum(df^2)+l)

df_pre=df

print(f(params))  

}

predict=c()

for(j in 1:(t-p+1)){  

y=c(1,Y[j:(j+p-1)]^2)

sig=sum(y*params)

predict=c(predict,sqrt(sig))

}

Ys=Y[(p):length(Y)]

plot(c(1:length(Ys)),Ys,xlim=c(1,length(Ys)),ylim=c(min(c(Ys,predict)),max(c(Ys,predict))),col=2,type="p",ylab="ARCH")

par(new=T)

plot(c(1:length(Ys)),predict,xlim=c(1,length(Ys)),ylim=c(min(c(Ys,predict)),max(c(Ys,predict))),col=3,type="p",ylab="ARCH")