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]]
______________________________________________
[email protected] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
