一般の時系列回帰過程

# AR多項式過程

library(dplyr)

data=data.frame(no=1:48,values=window(UKgas,c(1975,1),c(1999,3)))

r=4;dim=5

values=data$values

X=array(0,dim=c((nrow(data)-r),dim*r))

Y=c()

for(j in 1:(nrow(data)-r)){
  
vec=c()

for(i in 1:dim){
  
vec=c(vec,values[j:(j+r-1)]^i)  
  
}
  
X[j,]=vec

Y=c(Y,(values[r+j]))
  
}

X2=t(X)%*%X

sol_H2=eigen(X2)$vectors%*%diag(1/eigen(X2)$values)%*%t(eigen(X2)$vectors)

A=sol_H2%*%t(X)%*%Y

b=mean(Y)-sum(apply(X,2,mean)*A)

predict=X%*%A+b

result=data.frame(predict=predict,Y=Y)

N=nrow(X);

sigma_hat=sum(c(Y-X%*%A+b)^2)/N

AIC=log(sigma_hat)+2*r/N

# 罰則項付きのAR多項式過程

library(dplyr)

data=data.frame(no=1:48,values=window(UKgas,c(1975,1),c(1999,3)))

r=4;dim=5

values=data$values

X=array(0,dim=c((nrow(data)-r),dim*r))

Y=c()

for(j in 1:(nrow(data)-r)){
  
vec=c()

for(i in 1:dim){
  
vec=c(vec,values[j:(j+r-1)]^i)  
  
}
  
X[j,]=vec

Y=c(Y,(values[r+j]))
  
}

X2=t(X)%*%X


mesh=1

AIC_data=data.frame(lambda=seq(mesh,10000,mesh),AIC=0,CV=0,GCV=0)

for(i in 1:nrow(AIC_data)){

lambda=AIC_data$lambda[i]

sol_H=eigen(X2+diag(lambda,ncol(X)))$vectors%*%diag(1/eigen(X2+diag(lambda,ncol(X)))$values)%*%t(eigen(X2+diag(lambda,ncol(X)))$vectors)
  
A=sol_H%*%t(X)%*%Y

b=mean(Y)-sum(apply(X,2,mean)*A)

predict=X%*%A+b

sigma_hat=mean((Y-predict)^2)

H=X%*%sol_H%*%t(X)

AIC=length(Y)*log(2*pi+1)+length(Y)*log(sigma_hat)+2*sum(diag(H))

CV=mean(((Y-predict)/(1-diag(H)))^2)

GCV=mean(((Y-predict)/(1-mean(diag(H))))^2)

AIC_data$AIC[i]=AIC

AIC_data$GCV[i]=GCV

AIC_data$CV[i]=CV

}

plot(AIC_data$lambda,AIC_data$AIC)

lambda=min(AIC_data$lambda[AIC_data$AIC==min(AIC_data$AIC)])

sol_H=eigen(X2+diag(lambda,ncol(X)))$vectors%*%diag(1/eigen(X2+diag(lambda,ncol(X)))$values)%*%t(eigen(X2+diag(lambda,ncol(X)))$vectors)
  
A=sol_H%*%t(X)%*%Y

b=mean(Y)-sum(apply(X,2,mean)*A)

predict=X%*%A+b

result=data.frame(predict=predict,Y=Y)

# 一般の時系列回帰過程

library(dplyr)

data=data.frame(anscombe)

char=c("x1","x2","x3","x4")

char=c("x1")

r=4;cols=length(char)

values=data$y4

X=array(0,dim=c((nrow(data)-r),(1+cols)*r))


names=colnames(data)


Y=c()

for(j in 1:(nrow(data)-r)){

vec=c()

vec=c(vec,values[j:(j+r-1)])  

if(cols==1){

data_sub=array(data[,colnames(data) %in% char],dim=c(nrow(data),cols))

}else{
  
data_sub=data[,colnames(data) %in% char]  
  
}

for(i in 1:ncol(data_sub)){

vec=c(vec,data_sub[j:(j+r-1),i])  

}

X[j,]=vec

Y=c(Y,(values[r+j]))

}

X2=t(X)%*%X

sol_H2=eigen(X2)$vectors%*%diag(1/eigen(X2)$values)%*%t(eigen(X2)$vectors)

A=sol_H2%*%t(X)%*%Y

b=mean(Y)-sum(apply(X,2,mean)*A)

predict=X%*%A+b

result=data.frame(predict=predict,Y=Y)

N=nrow(X);

sigma_hat=sum(c(Y-X%*%A+b)^2)/N

AIC=log(sigma_hat)+2*r/N