[R] fitting functional autoregresive model

[Faheem Jan]

Hi, i trying to  extend the functional autoregressive model one FAR(1) to fit the functional autoregressive model of order two FAR(2). the coding i do for far(1) is library(fda)library(far)# CREATE DUMMY VARIABLESfactor2dummy=function(x){  n=length(x)  tab=as.factor(names(table(x)))  p=length(tab)  xdummy=matrix(0,n,p)  for(i in 1:p)  {    xdummy[x==tab[i],i]=1  }  colnames(xdummy)=tab  return(xdummy)}
# READ DATA
demnd=read.csv("c:/Users/Khan/Desktop/dem99141.csv",header=TRUE)xdata <- as.matrix(demnd[-1,-1], ncol = 25, nrow =1826, byrow= TRUE)class(xdata)date=strptime(as.character(xdata[,1]),"%Y-%m-%d")weekday=weekdays(date)week=format(date,"%U")xdata=xdata[,-1]xdata# DAILY AVERAGExmean=apply(xdata,1,mean)xmean
# SEASONAL ADJUSTMENT#seasadj=function(x) decompose(ts(x,frequency=7))$rand#xdata=apply(xdata,2,seasadj)#xdata=xdata[!is.na(xdata[,1]),]nrall=nrow(xdata)#wd=factor2dummy(weekday)#wnr=factor2dummy(week)#e=lm(xmean~wd+wnr-1)$residuals#tsdiag(arima(e,c(7,0,0)))#seasadj=function(x) lm(x~wd+wnr-1)$residuals#xdata=apply(xdata,2,seasadj)
# HOLD-OUT-PERIODnout=100nin=nrall-noutxin=xdata[1:nin,]
nr=nrow(xin)nc=ncol(xin)n=nr*ncy=matrix(t(xin),n,1)xfd=as.fdata(y,col=1,p=23,name="Cons") #p=23 is the multiple of length=39698 of data
# ESTIMATE FAR(1) MODELk1=far.cv(xfd,y="Cons",kn=20,ncv=1000)$minL2[1]far1=far(xfd,kn=k1)far1# ESTIMATE AR(1)-MODELSp=14f=function(x) ar(x,aic=FALSE,order.max=p)ar.models=apply(xin,2,f)
# FORECASTerrorsfar=matrix(0,nout,nc)errorsar=matrix(0,nout,nc)errorsnaive=matrix(0,nout,nc)predar=matrix(0,1,nc)prednaive=matrix(0,1,nc)for(i in 1:nout){  for(j in 1:length(ar.models))  {    predar[1,j]=predict(ar.models[[j]],newdata=xdata[(nr+i-p):(nr+i-1),j])$pred  }  xnew=as.fdata(t(xdata[nr+i-1,]),col=1,p=23,name="Cons")  pred=predict(far1,newdata=xnew)  prednaive=xdata[nr+i-1,]  obs=xdata[nr+i,]  errorsnaive[i,]=t(obs-prednaive)  errorsar[i,]=t(obs-predar)  errorsfar[i,]=t(obs-pred$Cons)}msefar=apply(errorsfar^2,2,mean)msefarmsear=apply(errorsar^2,2,mean)msenaive=apply(errorsnaive^2,2,mean)mean(msear)mean(msenaive)mean(msefar

i use the consumption data ...
	[[alternative HTML version deleted]]