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]]