sd_trim <- function(x,trim=0.2, const=TRUE){
# trimmed sd, where x is a matrix (column-wise)
x <- as.matrix(x)
if (const){
if (trim==0.1){const <- 0.7892}
else if (trim==0.2){const <- 0.6615}
else {warning("Did you specify the correct consistency constant for trimming?")}
}
else{const <- 1}
m <- apply(x,2,mean,trim)
res <- x-rep(1,nrow(x))%*%t(m)
qu <- apply(abs(res),2,quantile,1-trim)
sdtrim <- apply(matrix(res[t(abs(t(res))<=qu)]^2,ncol=ncol(x),byrow=FALSE),2,sum)
sdtrim <- sqrt(sdtrim/((nrow(x)*(1-trim)-1)))/const
return(sdtrim)
}
periodest<-function(x){
if (length(x)<20)
return (1)
x.spec <- spectrum(x,log="no",span=10,plot=FALSE)
spx <- x.spec$freq
spy <- 2*x.spec$spec
return(round(1/spx[which.max(spy)])) ## period estimation
}
nlcor<-function(x,y){
require(lazy)
N<-length(x)
I<-sample(1:N,round(N/3))
data<-data.frame(x[I],y[I])
y.lazy <- lazy(y ~ x,data,control=lazy.control(linIdPar=c(round(N/2),N)))
yh<-predict(y.lazy, newdata=x[-I])$h
cor(y[-I], yh)
}
#' timefit
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @description time fitting of time series
#' @details time fitting
#' @title timefit
#' @name timefit
#' @param TS.tr: training time series
#' @param C: number of fitted points
#' @param H: horizon
#' @return vector of H predictions
#' @export
timefit<-function(TS.tr,n,C,H){
Ntr=length(TS.tr)
Itr=seq(max(1,Ntr-(C+5)),Ntr)
Its=seq(Ntr+1,Ntr+H)
D=data.frame(cbind(Itr,TS.tr[Itr]))
names(D)<-c('t','ts')
weights=rev(exp(-(1:length(Itr))))
weights=weights/max(weights)
if (n==1)
mod=lm(ts~ 1,data=D,weights=weights)
if (n>1)
mod=lm(ts~ poly(t,min(length(Itr),n-1)),data=D,weights=weights)
Dts=data.frame(Its)
names(Dts)<-c('t')
TS.hat=predict(mod,newdata=Dts)
return(TS.hat)
}
#' KNN.multioutput
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \emph{Bontempi G. Ben Taieb S. Conditionally dependent strategies for multiple-step-ahead prediction in local learning, International Journal of Forecasting Volume 27, Issue 3, July–September 2011, Pages 689–699}
#' @references \url{https://tinyurl.com/sfmlh}
#' @description Multioutput KNN
#' @details Multioutput KNN which perform a locally constant model with weighted combination of local model on the basis of leave-one-out error. A tricube kernel is used to weight the distances.
#' @title KNN.multioutput
#' @name KNN.multioutput
#' @param X: training input [N,n]
#' @param Y: training output [N,m]
#' @param X.ts: test input [N.ts,n]
#' @param k: min number of neighbours
#' @param dist: type of distance: \code{euclidean, cosine}
#' @param F: forgetting factor
#' @param C: integer parameter which sets the maximum number of neighbours (Ck)
#' @param wta: if TRUE a winner-takes-all strategy is used; otherwise a weigthed combination is done on the basis of the l-o-o error
#' @param Reg: number (>1) of null terms to regularise the mean
#' @return vector of N.ts predictions
#' @export
#' @examples
#' ## Multi-step ahead time series forecasting
#' require(gbcode)
#' t=seq(0,200,by=0.1)
#' N<-length(t)
#' H<-50 ## horizon prediction
#' TS<-sin(t)
#' TS.tr=TS[1:(N-H)]
#' N.tr<-length(TS.tr)
#' TS.ts<-TS[(N-H+1):N]
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' E=MakeEmbedded(TS.tr,n=3,delay=0,hor=H,1)
#' X<-E$inp
#' Y<-E$out
#' N<-NROW(X)
#' Y.cont<-KNN.multioutput(X,Y,rev(TS.tr[(N.tr-H):N.tr]))
#' plot(t[(N-H+1):N],TS.ts)
#' lines(t[(N-H+1):N],Y.cont)
KNN.multioutput<- function(X,Y,X.ts,k=10,Di=NULL,
dist="euclidean",C=2,F=0,
wta=TRUE,scaleX=TRUE,
Reg=1){
Reg=max(Reg,1)
if (k<=0)
stop("k must be positive")
if (is.vector(X))
X<-array(X,c(length(X),1))
N<-NROW(X)
n<-NCOL(X)
if (is.vector(X.ts) & n>1){
N.ts<-1
X.ts<-array(X.ts,c(1,n))
} else {
if (n==1)
N.ts<-length(X.ts)
else
N.ts<-nrow(X.ts)
}
if (scaleX & all(apply(X,2,sd)>0.1)){
X=scale(X)
X.ts=scale(X.ts,attr(X,"scaled:center"), attr(X,"scaled:scale"))
}
if ( k >=NROW(X) ){
return (array(mean(Y),c(N.ts,1)))
}
if (is.vector(Y) ){
Y<-array(Y,c(length(Y),1))
}
m<-NCOL(Y)
if (n==1)
X<-array(X,c(N,1))
out.hat<-array(NA,c(N.ts,m))
if (is.null(Di)){
if (dist=="euclidean")
Ds<-dist2(X,X.ts)
if (dist=="cosine")
Ds<-dist.cos(X,X.ts)
} else
Ds<-Di
if (m>1){
w.na<-which(apply(Y,1,ana))
}else{
w.na<-which(is.na(Y))
}
if (length(w.na)>0)
Ds[w.na,]<-Inf
if (F <1)
Ds<-array(Ds*exp(-(F*(1:N))),c(N,N.ts)) ## forgetting factor
for (i in 1:N.ts){
index<-sort(Ds[,i],index.return=TRUE)
err<-numeric(C*k)+Inf
oo<-NULL
for (kk in k:(min(NROW(X),C*k))){
d<-Ds[index$ix[1:kk],i]/Ds[index$ix[kk+1],i]
## tricube kernel
wd<-c(((1-abs(d)^3)^3)*(abs(d)<1),numeric(Reg)+0.1)
wd<-wd/sum(wd)
if (any(is.na(wd)))
wd<-numeric(length(wd))+1/length(wd)
if (m>1){
#L<-apply(Y[index$ix[1:kk],],2,var)
L<-apply(Y[index$ix[1:kk],],2,constloo,wd[1:kk])
L<-L[which(is.finite(L))]
err[kk]<-mean(L)
YY=rbind(Y[index$ix[1:kk],],array(mean(Y),c(Reg,NCOL(Y))))
## add as many rows as Reg regularisation null terms
oo<-rbind(oo,apply(wd*YY,2,sum,na.rm=T))
} else {
err[kk]<-constloo(Y[index$ix[1:kk],1],wd[1:kk])
##err[kk]<-var(Y[index$ix[1:kk],1])
oo<-rbind(oo,mean(c(Y[index$ix[1:kk],1],numeric(Reg)+mean(Y)),na.rm=T))
}
if (is.na(err[kk])){
stop("KNN.multioutput error")
}
}
if (wta) {
if ((which.min(err)-k+1)<0)
stop("(which.min(err)-k+1)<0")
out.hat[i,]<-oo[which.min(err)-k+1,]
}else {
w<-(1/err[k:(C*k)])/sum(1/err[k:(C*k)])
## weights inversely proportional to loo errors
out.hat[i,]<-apply(oo*w,2,sum,na.rm=T)
}
}
out.hat
}
#' KNN.acf
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \emph{Bontempi G. Ben Taieb S. Conditionally dependent strategies for multiple-step-ahead prediction in local learning, International Journal of Forecasting Volume 27, Issue 3, July–September 2011, Pages 689–699}
#' @description Multioutput KNN
#' @details Multioutput KNN for multi-step-ahed prediction. It performs a locally constant model with weighted combination of local model on the basis of the dynamic properties of the training time series.
#' @title KNN.acf
#' @name KNN.acf
#' @param X: training input [N,n]
#' @param Y: training output [N,m]
#' @param X.ts: test input [N.ts,n]
#' @param k: min number of neighbours
#' @param dist: type of distance: \code{euclidean, cosine}
#' @param F: forgetting factor
#' @param C: integer parameter which sets the maximum number of neighbours (Ck)
#' @param wta: if TRUE a winner-takes-all strategy is used; otherwise a weigthed combination is done on the basis of the l-o-o error
#' @param Acf: autocorrelation function of the training series
#' @param Reg: number (>1) of null terms to regularise the mean
#' @return vector of N.ts predictions
#' @export
#' @examples
#' ## Multi-step ahead time series forecasting
#' require(gbcode)
#' library(lazy)
#' t=seq(0,200,by=0.1)
#' N<-length(t)
#' H<-500 ## horizon prediction
#' TS<-sin(t)+rnorm(N,sd=0.1)
#' TS.tr=TS[1:(N-H)]
#' N.tr<-length(TS.tr)
#' TS.ts<-TS[(N-H+1):N]
#' n=3
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' E=MakeEmbedded(TS.tr,n=n,delay=0,hor=H,1)
#' X<-E$inp
#' Y<-E$out
#' N<-NROW(X)
#' Y.cont<-KNN.acf(X,Y,rev(TS.tr[(N.tr-n+1):N.tr]),TS=TS.tr)
#' plot(t[(N-H+1):N],TS.ts)
#' lines(t[(N-H+1):N],Y.cont)
KNN.acf<- function(X,Y,X.ts,k=10,dist="euclidean",C=2,F=0,Acf,Pacf,TS,Reg=3){
Reg=max(Reg,1)
if (k<=0)
stop("k must be positive")
if (is.vector(X))
X<-array(X,c(length(X),1))
N<-NROW(X)
n<-NCOL(X)
if (is.vector(X.ts) & n>1){
N.ts<-1
X.tsts<-array(X.ts,c(1,n))
} else {
if ( n==1)
N.ts<-length(X.ts)
else
N.ts<-nrow(X.ts)
}
if ( k >=NROW(X)){
return (array(mean(Y),c(N.ts,1)))
}
m<-NCOL(Y)
if (n==1)
X<-array(X,c(N,1))
out.hat<-array(NA,c(N.ts,m))
Ds<-dist2(X,X.ts)
Ds[which(apply(Y,1,ana))]<-Inf
Ds<-array(Ds*exp(-(F*(1:N))),c(N,N.ts))
for (i in 1:N.ts){
index<-sort(Ds[,i],index.return=TRUE)
err<-numeric(C*k)+Inf
err2<-numeric(C*k)+Inf
oo<-NULL
aa<-NULL
ord.emb<-round(NCOL(X)/2)
ETS<-MakeEmbedded(array(TS,c(length(TS),1)),n=ord.emb,0)
for (kk in k:(C*k)){
d<-Ds[index$ix[1:kk],i]/Ds[index$ix[kk+1],i]
wd<-c(((1-abs(d)^3)^3)*(abs(d)<1),numeric(Reg)+0.1)
wd<-wd/sum(wd)
if (any(is.na(wd)))
wd<-1/length(wd)+numeric(length(wd))
if (m>1){
YY=rbind(Y[index$ix[1:kk],],array(mean(Y),c(Reg,NCOL(Y))))
LP<-apply(YY,2,mean,na.rm=T)
oo<-rbind(oo,LP)
ets<-MakeEmbedded(array(c(LP),c(length(LP),1)),n=ord.emb,0)
L<-apply(YY,2,constloo,wd)
L<-L[which(is.finite(L))]
err[kk]<-0
for (jj in 1:ord.emb) {
err[kk]<-err[kk]+mean((ets$out-
lazy.pred(X=ETS$inp[,1:jj],Y=ETS$out,
X.ts=rbind(ets$inp[,1:jj]),linPar=c(5,10),class=FALSE))^2)
}
err2[kk]<-mean(L)
} else {
stop("Only for multi step ahead prediction")
}
if (is.na(err[kk]))
stop("is.na(err[kk])")
}
w2<-(1/err2[k:(C*k)])/sum(1/err2[k:(C*k)])
out.hat[i,]<-oo[which.min(err)-k+1,]
}
out.hat
}
#' KNN.acf.lin
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \emph{Bontempi G. Ben Taieb S. Conditionally dependent strategies for multiple-step-ahead prediction in local learning, International Journal of Forecasting Volume 27, Issue 3, July–September 2011, Pages 689–699}
#' @description Multioutput KNN
#' @details Multioutput KNN for multi-step-ahed prediction. It performs a locally constant model with weighted combination of local model on the basis of the autocorrelation and partial correlation properties of the training time series.
#' @title KNN.acf.lin
#' @name KNN.acf.lin
#' @param X: training input [N,n]
#' @param Y: training output [N,m]
#' @param X.ts: test input [N.ts,n]
#' @param k: min number of neighbours
#' @param dist: type of distance: \code{euclidean, cosine}
#' @param F: forgetting factor
#' @param C: integer parameter which sets the maximum number of neighbours (Ck)
#' @param wta: if TRUE a winner-takes-all strategy is used; otherwise a weigthed combination is done on the basis of the l-o-o error
#' @param Acf: autocorrelation function of the training series
#' @param Reg: number (>1) of null terms to regularise the mean
#' @return vector of N.ts predictions
#' @export
#' @examples
#' ## Multi-step ahead time series forecasting
#' require(gbcode)
#' t=seq(0,200,by=0.1)
#' N<-length(t)
#' H<-500 ## horizon prediction
#' TS<-sin(t)+rnorm(N,sd=0.1)
#' TS.tr=TS[1:(N-H)]
#' N.tr<-length(TS.tr)
#' TS.ts<-TS[(N-H+1):N]
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' E=MakeEmbedded(TS.tr,n=3,delay=0,hor=H,1)
#' X<-E$inp
#' Y<-E$out
#' N<-NROW(X)
#' ACF.lag<-5
#' Y.cont<-KNN.acf.lin(X,Y,rev(TS.tr[(N.tr-H):N.tr]),Acf=acf(TS.tr,lag.max=ACF.lag,plot=FALSE)$acf,Pacf=pacf(TS.tr,lag.max=ACF.lag,plot=FALSE)$acf,TS=TS.tr)
#' plot(t[(N-H+1):N],TS.ts)
#' lines(t[(N-H+1):N],Y.cont)
KNN.acf.lin<- function(X,Y,X.ts,k=10,dist="euclidean",C=2,F=0,Acf,Pacf,TS,Reg=3){
Reg=max(Reg,1)
if (k<=0)
stop("k must be positive")
if (is.vector(X))
X<-array(X,c(length(X),1))
N<-NROW(X)
n<-NCOL(X)
if (is.vector(X.ts) & n>1){
N.ts<-1
X.ts<-array(X.ts,c(1,n))
} else {
if (n==1)
N.ts<-length(X.ts)
else
N.ts<-nrow(X.ts)
}
if ( k >=NROW(X)){
return (array(mean(Y),c(N.ts,1)))
}
m<-NCOL(Y)
if (n==1)
X<-array(X,c(N,1))
out.hat<-array(NA,c(N.ts,m))
Ds<-dist2(X,X.ts)
Ds[which(apply(Y,1,ana))]<-Inf
Ds<-array(Ds*exp(-(F*(1:N))),c(N,1))
for (i in 1:N.ts){
index<-sort(Ds[,i],index.return=TRUE)
err<-numeric(C*k)+Inf
err2<-numeric(C*k)+Inf
oo<-NULL
for (kk in k:(C*k)){
d<-Ds[index$ix[1:kk],i]/Ds[index$ix[kk+1],i]
##d<-numeric(length(d))
##wd<-exp(-d^2)
wd<-c(((1-abs(d)^3)^3)*(abs(d)<1),numeric(Reg)+0.1)
wd<-wd/sum(wd)
if (any(is.na(wd)))
wd<-1/length(wd)+numeric(length(wd))
if (m>1){
YY=rbind(Y[index$ix[1:kk],],array(mean(Y),c(Reg,NCOL(Y))))
LP<-apply(YY,2,mean,na.rm=T)
oo<-rbind(oo,LP)
A<-acf(c(TS,LP),lag.max=length(Acf)-1,plot=FALSE)$acf
PA<-pacf(c(TS,LP),lag.max=length(Pacf),plot=FALSE)$acf
err[kk]<- 1-abs(cor(c(PA),c(Pacf)))+1-abs(cor(c(A),c(Acf)))
} else {
O<-mean(c(Y[index$ix[1:kk],1],numeric(Reg)+mean(Y)),na.rm=T)
A<-acf(c(O),lag.max=length(Acf)-1,plot=FALSE)$acf
PA<-pacf(c(O),lag.max=length(Pacf),plot=FALSE)$acf
err[kk]<-1-abs(cor(c(PA),c(Pacf)))##+1-abs(cor(c(A),c(Acf)))
wwP<-exp(-0.5*(1:length(PA)))
wwA<-exp(-0.5*(1:length(A)))
err[kk]<-mean((wwP*(c(PA)-c(Pacf)))^2)+mean(wwA*((c(A)-c(Acf)))^2)
oo<-rbind(oo,0)
}
if (is.na(err[kk]))
stop("Error")
}
out.hat[i,]<-oo[which.min(err)-k+1,]
}
out.hat
}
#' KNN.pls
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \emph{Bontempi G. Ben Taieb S. Conditionally dependent strategies for multiple-step-ahead prediction in local learning, International Journal of Forecasting Volume 27, Issue 3, July–September 2011, Pages 689–699}
#' @description Multioutput KNN
#' @details Multioutput KNN for multi-step-ahed prediction. It performs a locally constant model with weighted combination of local model on the basis of partial least-squares error
#' @title KNN.pls
#' @name KNN.pls
#' @param X: training input [N,n]
#' @param Y: training output [N,m]
#' @param X.ts: test input [N.ts,n]
#' @param k: min number of neighbours
#' @param dist: type of distance: \code{euclidean, cosine}
#' @param F: forgetting factor
#' @param C: integer parameter which sets the maximum number of neighbours (Ck)
#' @param wta: if TRUE a winner-takes-all strategy is used; otherwise a weigthed combination is done on the basis of the l-o-o error
#' @param Acf: autocorrelation function of the training series
#' @return vector of N.ts predictions
#' @export
#' @examples
#' ## Multi-step ahead time series forecasting
#' library(pls)
#' require(gbcode)
#' t=seq(0,400,by=0.1)
#' N<-length(t)
#' H<-1500 ## horizon prediction
#' TS<-sin(t)+rnorm(N,sd=0.1)
#' TS.tr=TS[1:(N-H)]
#' N.tr<-length(TS.tr)
#' TS.ts<-TS[(N-H+1):N]
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' n=3
#' E=MakeEmbedded(TS.tr,n=n,delay=0,hor=H,1)
#' X<-E$inp
#' Y<-E$out
#' ACF.lag<-5
#' Y.cont<-KNN.pls(X,Y,rev(TS.tr[(N.tr-n+1):N.tr]))
#' plot(t[(N-H+1):N],TS.ts)
#' lines(t[(N-H+1):N],Y.cont,col="red")
#'
#'
#' rm(list=ls())
#' require(gbcode)
#' data(A)
#' N<-NROW(A)
#' H<-200 ## horizon prediction
#' TS<-A[,1]
#' TS.tr=TS[1:(N-H)]
#' N.tr<-length(TS.tr)
#' TS.ts<-TS[(N-H+1):N]
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' n=16
#' E=MakeEmbedded(TS.tr,n=16,delay=0,hor=H,1)
#' X<-E$inp
#' Y<-E$out
#' Y.cont<-KNN.pls(X,Y,rev(TS.tr[(N.tr-n+1):N.tr]))
#' plot((N-H+1):N,TS.ts)
#' lines((N-H+1):N,Y.cont)
#'
KNN.pls<- function(X,Y,X.ts,k=10,dist="euclidean",C=2,F=0,del=0){
if (k<=0)
stop("k must be positive")
if (is.vector(X))
X<-array(X,c(length(X),1))
N<-NROW(X)
n<-NCOL(X)
if (is.vector(X.ts) & n>1){
N.ts<-1
X.ts<-array(X.ts,c(1,n))
} else {
if (n==1)
N.ts<-length(X.ts)
else
N.ts<-nrow(X.ts)
}
if ( k >=NROW(X)){
return (array(mean(Y),c(N.ts,1)))
}
m<-NCOL(Y)
if (n==1)
X<-array(X,c(N,1))
out.hat<-array(NA,c(N.ts,m))
init<-1
Ds<-dist2(X[,init:n],X.ts[,init:n])
Ds[which(apply(Y,1,ana))]<-Inf
Ds<-array(Ds*exp(-(F*(1:N))),c(N,1))
for (i in 1:N.ts){
index<-sort(Ds[,i],index.return=TRUE)
err<-numeric(C*k)+Inf
err2<-numeric(C*k)+Inf
oo<-NULL
for (kk in k:(C*k)){
d<-Ds[index$ix[1:kk],i]/Ds[index$ix[kk+1],i]
XX<-data.frame(X[index$ix[1:kk],])
names(XX)<-as.character(1:NCOL(XX))
XXTs<-data.frame(array(X.ts[i,],c(1,n)))
names(XXTs)<-as.character(1:NCOL(XX))
YY<-Y[index$ix[1:kk],]
MV<-plsr(YY~.,data=XX,validation="LOO")
LP<-predict(MV,newdata=XXTs)
nc<-which.min(apply(MV$validation$PRESS,2,mean,na.rm=T))
LP<-c(LP[,,nc])
oo<-rbind(oo,LP)
err[kk]<- mean(MV$validation$PRESS[,nc],na.rm=T)
if (is.na(err[kk]))
stop("Error")
}
w2<-(1/err[k:(C*k)])/sum(1/err[k:(C*k)])
out.hat[i,]<-apply(oo*w2,2,sum,na.rm=T)
}
out.hat
}
#' lin.pls
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \emph{Bontempi G. Ben Taieb S. Conditionally dependent strategies for multiple-step-ahead prediction in local learning, International Journal of Forecasting Volume 27, Issue 3, July–September 2011, Pages 689–699}
#' @description Multioutput lin PLS
#' @details Multioutput PLS for multi-step-ahed prediction. It performs a partial least-squares
#' @title lin.pls
#' @name lin.pls
lin.pls<- function(X,Y,X.ts){
if (is.vector(X))
X<-array(X,c(length(X),1))
N<-NROW(X)
n<-NCOL(X)
if (is.vector(X.ts) & n>1){
N.ts<-1
X.ts<-array(X.ts,c(1,n))
} else {
if (n==1)
N.ts<-length(X.ts)
else
N.ts<-nrow(X.ts)
}
m<-NCOL(Y)
if (n==1)
X<-array(X,c(N,1))
out.hat<-array(NA,c(N.ts,m))
XX<-data.frame(X)
names(XX)<-as.character(1:NCOL(XX))
XXTs<-data.frame(X.ts)
names(XXTs)<-as.character(1:NCOL(XX))
MV<-plsr(Y~.,data=XX,validation="CV")
LP<-predict(MV,newdata=XXTs)
nc<-which.min(apply(MV$validation$PRESS,2,mean,na.rm=T))
out.hat<-c(LP[,,nc])
}
#' multiplestepAhead
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \emph{Bontempi G. Ben Taieb S. Conditionally dependent strategies for multiple-step-ahead prediction in local learning, International Journal of Forecasting Volume 27, Issue 3, July–September 2011, Pages 689–699}
#' @references \url{https://tinyurl.com/sfmlh}
#' @description multiplestepAhead: : univariate multi-step-ahead forecaster
#' @details Wrapper over a set of methods for univariate multi-step-ahead time series forecasting
#' @title multiplestepAhead
#' @name multiplestepAhead
#' @param TS: time series in column format [T,1]
#' @param n: embedding order
#' @param H: horizon
#' @param dist: type of distance: \code{euclidean, cosine} for lazy methods
#' @param F: forgetting factor
#' @param C: integer parameter which sets the maximum number of neighbours (C*k) for lazy methods
#' @param detrend: real parameter (in [0,1]) which fixes the size of window used for linear detrending the series (0 corresponds to all series and 1 to ten latest terms). If detrend<0 no detrending is carried out
#' @param smooth: if TRUE, the prediction is obtained by averaging multiple windows with different starting points
#' @param engin: if TRUE, a number of additional features (related to the quantiles) are engineered and added
#' @param method:
#' \itemize{
#' \item{arima}: prediction based on the \pkg{forecast} package
#' \item{stat_naive}: naive predictor based on the M4 competition code
#' \item{stat_ses_naive}: prediction based on the M4 competition code
#' \item{stat_naive2}: naive predictor based on the M4 competition code
#' \item{stat_ses}: SES predictor based on the M4 competition code
#' \item{stat_holt}: Holt predictor based on the M4 competition code
#' \item{stat_damped}: prediction based on the M4 competition code
#' \item{stat_theta}: Theta predictor based on the M4 competition code
#' \item{stat_comb}: prediction based on the M4 competition code
#' \item{direct}: direct prediction based on \link{KNN.multioutput} function
#' \item{iter}: recursive prediction based on \link{KNN.multioutput} function
#' \item{lazydirect}: locally linear direct prediction based on \link{lazy.pred} function
#' \item{clazydirect}: locally constant direct prediction based on \link{lazy.pred} function
#' \item{lazyiter}: recursive prediction based on \link{lazy.pred} function
#' \item{lindirect}: direct prediction based on linear predictor
#' \item{rfdirect}: direct prediction based on Random Forest predictor
#' \item{rfiter}: recursive prediction based on Random Forest predictor
#' \item{mimo}: MIMO prediction based on \link{KNN.multioutput} function
#' \item{mimo.comb}: MIMO prediction based on \link{KNN.multioutput} function which combines a set of predictors based on different horizons and different starting points
#' \item{mimo.acf}: MIMO prediction based on \link{KNN.acf} function which combines a set of predictors based on different horizons and different starting points
#' \item{mimo.acf.lin}: MIMO prediction based on \link{KNN.acf.lin} function which combines a set of predictors based on different horizons and different starting points
#' \item{mimo.pls}: MIMO prediction based on \link{KNN.pls} function which combines a set of predictors based on different horizons and different starting points
#' \item{mimo.lin.pls}: MIMO prediction based on Partial Least Squares which combines a set of predictors based on different horizons and different starting points
#' \item{mimo_rr}: MIMO prediction based on linear ridge regression
#' \item{mimo_red}: MIMO prediction based on linear reduced rank regression (\pkg{rrpack})
#' \item{mimo_las}: MIMO prediction based on linear lasso regression based on python (\pkg{reticulate})
#' \item{mimo_cca}: MIMO prediction based on linear canonical correlation
#' \item{mimo_red}: MIMO prediction based on linear reduced rank regression
#' \item{rnn}: prediction based on python (\pkg{reticulate}) implementation of rnn (recurrent neural networks)
#' \item{lstm}: prediction based on python (\pkg{reticulate}) implementation of lstm neural networks
#' \item{transf}: prediction based on python (\pkg{reticulate}) implementation of transformer neural networks
#' }
#' @return H step ahead predictions
#' @details
#' The python forecasters require the installation of \pkg{reticulate} and several python packages (scikit-learn, tensorflow, keras)
#' @export
#' @examples
#' ## Multi-step ahead time series forecasting
#'
#' rm(list=ls())
#' t=seq(0,400,by=0.1)
#' N<-length(t)
#' H<-500 ## horizon prediction
#' TS<-sin(t)+rnorm(N,sd=0.1)
#' TS.tr=TS[1:(N-H)]
#' N.tr<-length(TS.tr)
#' TS.ts<-TS[(N-H+1):N]
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' Y.cont=multiplestepAhead(TS.tr,n=3, H=H,method="mimo")
#'
#' plot(t[(N-H+1):N],TS.ts)
#' lines(t[(N-H+1):N],Y.cont,col="red")
#'
#'
#'
#'
#' ## Multi-step ahead time series forecasting Santa Fe chaotic time series A
#' rm(list=ls())
#' require(gbcode)
#' data(A)
#' TS=A
#' N<-1000
#' H<-200
#' TS.tr=TS[1:N,1]
#' TS.ts<-TS[(N+1):(N+H),1]
#' TS.tr=array(TS.tr,c(length(TS.tr),1))
#' Y.dir=multiplestepAhead(TS.tr,n=12, H=H,method="lazydirect")
#' Y.mimo.comb=multiplestepAhead(TS.tr,n=12, H=H,method="mimo.comb")
#' plot((N-H+1):N,TS.ts,type="l",xlab="",ylab="Santa Fe A series")
#' lines((N-H+1):N,Y.dir,col="red")
#' lines((N-H+1):N,Y.mimo.comb,col="green")
#'
#'
#'
multiplestepAhead<-function(TS,n,H,D=0, method="direct",
FF=0,smooth=FALSE,maxfs=6,
XC=NULL,detrend=0, forget=-1, engin=FALSE,
Kmin=5,C=2,debug=FALSE, learner="rf",
verbose=FALSE,...){
args<-list(...)
if (length(args)>0)
for(i in 1:length(args)) {
assign(x = names(args)[i], value = args[[i]])
}
if (NCOL(TS)>1)
stop("Only for univariate time series")
if (any(is.na(TS)))
stop("NA in the series")
N<-length(TS)
if (H>round(N/2))
stop("Too long horizon ")
qu=seq(0.1,1,by=0.1)
if ((N-n-H)<10)
method="stat_naive"
if (forget>0){
I=min(N-10,max(1,round(N*forget))):N
TS=TS[I]
if (!is.null(XC))
XC=XC[I[1]:NROW(XC),]
N<-length(TS)
}
trnd.ts=numeric(H)
trnd.ts2=numeric(H)
if (detrend>0){
S=detectSeason(TS,Ls=N+H,pmin=detrend)
TS=TS-S$spattern[1:N]-S$strend[1:N]
trnd.ts<-S$spattern+S$strend
trnd.ts<-trnd.ts[(N+1):(N+H)]
}
if (!is.null(XC)){
trnd2=pred("lin",XC[1:N,],TS,XC[1:N,],classi=FALSE,lambda=1e-2)
trnd.ts2=pred("lin",XC[1:N,],TS,XC,classi=FALSE,lambda=1e-2)[(N+1):(N+H)]
TS=TS-trnd2
trnd.ts=trnd.ts+trnd.ts2
}
### Set of statistical methods borrowed from FPP3
if (method=="fpp_ETS"){
library(fpp3)
N=length(TS)
TSI=tsibble(
date = as.Date(Sys.Date()) + 1:N,
index='date',
value = TS
)
fit<-TSI%>%model(ETS(value))
p=data.frame(fit%>%forecast(h=H))[,".mean"]
return(c(p+trnd.ts))
}
if (method=="fpp_Holt"){
library(fpp3)
N=length(TS)
TSI=tsibble(
date = as.Date(Sys.Date()) + 1:N,
index='date',
value = TS
)
fit<-TSI%>%model(model1=ETS(value~error("A")+trend("A")+season("N")))
p=data.frame(fit%>%forecast(h=H))[,".mean"]
return(c(p+trnd.ts))
}
if (method=="fpp_Damped"){
library(fpp3)
N=length(TS)
TSI=tsibble(
date = as.Date(Sys.Date()) + 1:N,
index='date',
value = TS
)
fit<-TSI%>%model(model1=ETS(value~error("A")+trend("Ad")+season("N")))
p=data.frame(fit%>%forecast(h=H))[,".mean"]
return(c(p+trnd.ts))
}
if (method=="fpp_ARIMA"){
library(fpp3)
N=length(TS)
TSI=tsibble(
date = as.Date(Sys.Date()) + 1:N,
index='date',
value = TS
)
fit<-TSI%>%model(model1=ARIMA(value))
p=data.frame(fit%>%forecast(h=H))[,".mean"]
return(c(p+trnd.ts))
}
if (method=="fpp_prophet"){
library(fpp3)
library(fable.prophet)
N=length(TS)
TSI=tsibble(
date = as.Date(Sys.Date()) + 1:N,
index='date',
value = TS
)
fit<-TSI%>%model(model1=prophet(value))
p=data.frame(fit%>%forecast(h=H))[,".mean"]
return(c(p+trnd.ts))
}
if (method=="fpp_NNETAR"){
library(fpp3)
N=length(TS)
TSI=tsibble(
date = as.Date(Sys.Date()) + 1:N,
index='date',
value = TS
)
fit<-TSI%>%model(model1=NNETAR(value))
p=data.frame(fit%>%forecast(h=H))[,".mean"]
return(c(p+trnd.ts))
}
### Set of statistical methods borrowed from M4 competition
if (method=="stat_naive"){
p=StatPredictors1(TS, H , index=1,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_ses_naive"){
p=StatPredictors1(TS, H , index=2,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_naive2"){
p=StatPredictors1(TS, H , index=3,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_ses"){
p=StatPredictors1(TS, H , index=4,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_holt"){
p=StatPredictors1(c(TS), H , index=5,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_damped"){
p=StatPredictors1(c(TS), H , index=6,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_theta"){
p=StatPredictors1(c(TS), H , index=7,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_avg"){
p=numeric(H)+mean(c(TS))
return(c(p+trnd.ts))
}
if (method=="stat_comb"){
p=StatPredictors1(c(TS), H , index=8,verbose=F)
return(c(p+trnd.ts))
}
if (method=="stat_4theta"){
p=StatPredictors1(c(TS), H , index=9,verbose=F)
return(c(p+trnd.ts))
}
if (method=="mimo_rr"){
p=multiridge(array(TS,c(length(TS),1)),n,H,MIMO=TRUE,verbose=verbose,...)$Yhat
return(c(p+trnd.ts))
}
if (method=="mimo_red"){
p=multirr(array(TS,c(length(TS),1)),n,H,B=0,QRdec=FALSE,verbose=verbose,...)
return(c(p+trnd.ts))
}
if (method=="mimo_cca"){
p=multicca(array(TS,c(length(TS),1)),n,H,...)
return(c(p+trnd.ts))
}
if (method=="mimo_las"){
p=multiml(cbind(TS),n,H,learner="py.lasso_regr",...)
return(c(p+trnd.ts))
}
### keras based RNN: it requires keras
if (method=="rnn"){
p=pyrnnpredgpt(cbind(TS),H,n=n,...)
##p=rnnpred(array(TS,c(length(TS),1)), H,...)
return(c(p+trnd.ts))
}
### keras based RNN: it requires keras
if (method=="lstm"){
p=pylstmpredgpt(cbind(TS),H,n=n,...)
#p=lstmpred(array(TS,c(length(TS),1)), H,...)
return(c(p+trnd.ts))
}
if (method=="transf"){
p=pytransfpredgpt(cbind(TS),H,n=n,...)
return(c(p+trnd.ts))
}
if (method=="nbeats"){
p=pytorchnbeats(cbind(TS),H,n=n,...)
return(c(p+trnd.ts))
}
if (method=="deepar"){
p=pytorchndeepar(cbind(TS),H,n=n,...)
return(c(p+trnd.ts))
}
if (method=="tft"){
p=pytorchntft(cbind(TS),H,n=n,...)
return(c(p+trnd.ts))
}
if (method=="gluon"){
p=gluonpred(array(TS,c(length(TS),1)),n=n, H=H,...)
return(c(p+trnd.ts))
}
if (method=="nbeats"){
p=nbeatspred(array(TS,c(length(TS),1)),n=n, H=H,...)
return(c(p+trnd.ts))
}
if (method=="darts_nbeats"){
p=dartsnbeats(array(TS,c(length(TS),1)),n=n, H=H,...)
return(c(p+trnd.ts))
}
if (method=="darts_tft"){
p=dartstft(array(TS,c(length(TS),1)),n=n, H=H,...)
return(c(p+trnd.ts))
}
if (method=="nbeatsens"){
p=nbeatsenspred(array(TS,c(length(TS),1)),n=n, H=H,...)
return(c(p+trnd.ts))
}
if (sd_trim(TS)<0.001 && method !="timefit" )
return (numeric(H)+TS[1]+trnd.ts)
TS<-array(TS,c(N,1))
M <- MakeEmbedded(TS, n, D, H, w = 1)
X<-M$inp
Y<-M$out
if (n==1)
X<-array(X,c(length(X),1))
NX=NROW(X)
select.var=1:NCOL(X)
if (length(select.var)>maxfs ) {
if (engin){
X<-cbind(X,t(apply(X,1,quantile,qu)))
}
rfs=numeric(NCOL(X))
for (j in 1:NCOL(Y)){
fs=mrmr(X,Y[,j],nmax=min(NCOL(X)-1,maxfs))
rfs[fs]=rfs[fs]+1
}
select.var=sort(rfs,decr=TRUE,index=TRUE)$ix[1:min(maxfs,NCOL(X)-1)]
}
if (length(select.var)>2){
sdX<-apply(X[,select.var],2,sd)
ws<-which(sdX<0.01)
if (length(ws)>0)
select.var<-setdiff(select.var,ws)
}
if (length(select.var)<1)
return(numeric(H)+mean(Y))
if (length(select.var)==1)
if (sd(X)<0.01)
return(numeric(H)+mean(Y))
q<-TS[seq(N-D,N-n+1-D,by=-1),1]
if (engin)
q=c(q,quantile(q,qu))
## TS=[TS(1), TS(2),....., TS(N)]
## D=0: q=[TS(N), TS(N-1),...,TS(N-n+1)]
switch(method,
arima={
fit <- tryCatch(
{
arima(TS,arimaorder(auto.arima(TS)))
},
error = function(e){
arima(TS,c(0,D,1))
}
)
p<-predict(fit,H)$pred
},
timefit={
p<-timefit(TS,n,C,H)
},
direct={
p<-numeric(H)
for (h in 1:H){
I<-1:(NROW(X))
Yh=Y[,h]
if (length(which(!is.na(Yh)))<1)
p[h]=0
else {
if (length(which(!is.na(Yh)))<NROW(X) ){
p[h]<-mean(Yh,na.rm=TRUE)
}else{
p[h]<-KNN.multioutput(X[,select.var],array(Yh,c(NX,1)),
q[select.var],k=Kmin,C=C,F=FF,Reg=1)
}
}
}
},
lazydirect={
p<-numeric(H)
for (h in 1:H){
wna=which(!is.na(Y[,h]))
if (length(wna)<1){
p[h]=0
}else{
if (length(wna)>9){
Xw=array(X[wna,select.var],c(length(wna),length(select.var)))
Yw=array(Y[wna,h],c(length(wna),1))
LPar=c(Kmin,(C+1)*Kmin)*length(select.var)
LPar[1]=min(LPar[1],NROW(Xw)-2)
LPar[2]=min(LPar[2],NROW(Xw)-1)
CPar=c(Kmin,C*Kmin+1)
CPar[1]=min(CPar[1],NROW(Xw)-1)
CPar[2]=min(CPar[2],NROW(Xw))
if (NROW(Xw) <= (7*NCOL(Xw)))
LPar=NULL
CPar=c(Kmin,C*Kmin)
CPar[1]=min(CPar[1],NROW(X)-1)
vX=1
if (!is.vector(Xw))
vX=apply(Xw,2,sd)
if (any(vX<0.01)){
v0=which(vX<0.01)
if (length(v0)>= length(select.var)){
p[h]=mean(Y[,h],na.rm=TRUE)
} else {
q2=q[select.var]
p[h]<-lazy.pred(Xw[,-v0],Yw,rbind(q2[-v0]),
conPar=CPar,linPar=LPar,cmbPar=3)
}
}else{
p[h]<-lazy.pred(Xw,Yw,rbind(q[select.var]),
conPar=CPar,linPar=LPar,cmbPar=3)
}
}else
p[h]=mean(Y[,h],na.rm=TRUE)
}
} ## for h
},
clazydirect={
p<-numeric(H)
for (h in 1:H){
wna=which(!is.na(Y[,h]))
if (length(wna)<1){
p[h]=0
}else{
if (length(wna)>9){
Xw=X[wna,select.var]
Yw=array(Y[wna,h],c(length(wna),1))
CPar=c(Kmin,C*Kmin+1)
CPar[1]=min(CPar[1],NROW(Xw)-1)
CPar[2]=min(CPar[2],NROW(Xw))
LPar=NULL
vX=1
if (!is.vector(Xw))
vX=apply(Xw,2,sd)
if (any(vX<0.01)){
v0=which(vX<0.01)
if (length(v0)>= length(select.var)){
p[h]=mean(Y[,h],na.rm=TRUE)
} else {
q2=q[select.var]
p[h]<-lazy.pred(Xw[,-v0],Yw,rbind(q2[-v0]),
conPar=CPar,linPar=LPar,cmbPar=3)
}
}else{
p[h]<-lazy.pred(Xw,Yw,rbind(q[select.var]),
conPar=CPar,linPar=LPar,cmbPar=3)
}
}else
p[h]=mean(Y[,h],na.rm=TRUE)
}
} ## for h
},
rfdirect={
p<-numeric(H)
for (h in 1:H){
wna=which(!is.na(Y[,h]))
if (length(wna)<1){
p[h]=0
}else{
if (length(wna)>9){
p[h]<-rf.pred(X[wna,select.var],array(Y[wna,h],c(length(wna),1)),q[select.var],
class=FALSE,ntree=C*50)
}else
p[h]=mean(Y[,h],na.rm=TRUE)
}
}
},
mldirect={
p<-numeric(H)
for (h in 1:H){
wna=which(!is.na(Y[,h]))
if (length(wna)<1){
p[h]=0
}else{
if (length(wna)>9){
p[h]<-pred(learner,X[wna,select.var],array(Y[wna,h],c(length(wna),1)),q[select.var],
class=FALSE,ntree=C*50)
}else
p[h]=mean(Y[,h],na.rm=TRUE)
}
}
},
lindirect={
p<-numeric(H)
for (h in 1:H){
wna=which(!is.na(Y[,h]))
if (length(wna)<1){
p[h]=0
}else{
if (length(wna)>9){
p[h]<-lin.pred(X[wna,select.var],array(Y[wna,h],c(length(wna),1)),q[select.var],
class=FALSE,lambda=1e-3*C)
}else
p[h]=mean(Y[,h],na.rm=TRUE)
}
} ## for h
},
mimo={
p<-KNN.multioutput(X[,select.var],Y,q[select.var],k=Kmin,C=C,F=FF,Reg=1)
},
mimo.comb={
pdirect2<-NULL
if (smooth & H >=2) ## start before the end of the series with an horizon H
for (h in round(H/2):(H)){
p2<-numeric(H)+NA
q2<-TS[seq(N-H+h-D,N+1-n-H+h-D,by=-1),1]
if (engin)
q2=c(q2,quantile(q2,qu))
## TS=[TS(1), TS(2),....., TS(N)]
## D=0: q2=[TS(N-H+h), TS(N-1-H+h),...,TS(N-n+1-H+h)]
## pred= [TS(N-H+h+1),...TS(N+h+1)]
KK<-KNN.multioutput(X[,select.var],Y,q2[select.var],k=Kmin,C=C,F=FF)
p2[1:h]<-KK[(H-h+1):H]
pdirect2<-rbind(pdirect2,p2)
}
for (h in round(H/2):(H)){ ## start at the end of the series with different horizons h=[H/2,...H]
p2<-numeric(H)+NA
q2<-TS[seq(N-D,N+1-n-D,by=-1),1]
if (engin)
q2=c(q2,quantile(q2,qu))
KK<-KNN.multioutput(X[,select.var],Y[,1:h],q2[select.var],k=Kmin,C=C,F=FF)
p2[1:h]<-KK
pdirect2<-rbind(pdirect2,p2)
}
## combination of different predictions
p<-apply(pdirect2,2,mean,na.rm=T)
},
mimo.acf={
pdirect3<-NULL
TS.acf<-TS[,1]
if (smooth & H >=2)
for (h in round(H/2):(H)){ ## start before the end of the series with an horizon H
p2<-numeric(H)+NA
q2<-TS[seq(N-H+h-D,N+1-n-H+h-D,by=-1),1]
KK<-KNN.acf(X[,select.var],Y,q2[select.var],k=Kmin,C=C,F=FF,
TS=TS.acf,D,Reg=1)
p2[1:h]<-KK[(H-h+1):H]
pdirect3<-rbind(pdirect3,p2)
}
for (h in round(H/2):(H)){
p2<-numeric(H)+NA
q2<-TS[seq(N-D,N+1-n-D,by=-1),1]
KK<-KNN.acf(X[,select.var],Y[,1:h],q2[select.var],k=Kmin,C=C,F=FF,
TS=TS.acf,D,Reg=1)
p2[1:h]<-KK
pdirect3<-rbind(pdirect3,p2)
}
## combination of different predictions
p<-apply(pdirect3,2,mean,na.rm=T)
},
mimo.acf.lin={
ACF.lag<-5
pdirect4<-NULL
TS.acf<-TS ##i-1=N
if (smooth & H >=2)
for (h in round(H/2):(H)){ ## start before the end of the series with an horizon H
p2<-numeric(H)+NA
q2<-TS[seq(N-H+h-D,N+1-n-H+h-D,by=-1),1]
KK<-KNN.acf.lin(X[,select.var],Y,q2[select.var],k=Kmin,C=C,F=FF,
Acf=acf(TS.acf,lag.max=ACF.lag,plot=F)$acf,
Pacf=pacf(TS.acf,lag.max=ACF.lag,plot=F)$acf,TS=TS.acf,D,
Reg=1)
p2[1:h]<-KK[(H-h+1):H]
pdirect4<-rbind(pdirect4,p2)
}
for (h in round(H/2):(H)){
p2<-numeric(H)+NA
q2<-TS[seq(N-D,N+1-n-D,by=-1),1]
KK<-KNN.acf.lin(X[,select.var],Y[,1:h],q2[select.var],k=Kmin,C=C,F=FF,
Acf=acf(TS.acf,lag.max=ACF.lag,plot=F)$acf,
Pacf=pacf(TS.acf,lag.max=ACF.lag,plot=F)$acf,TS=TS.acf,D,Reg=1)
p2[1:h]<-KK
pdirect4<-rbind(pdirect4,p2)
}
## combination of different predictions
p<-apply(pdirect4,2,mean,na.rm=T)
},
mimo.pls={
pdirect5<-NULL
if (smooth & H >=2)
for (h in round(H/2):(H)){ ## start before the end of the series with an horizon H
p2<-numeric(H)+NA
q2<-TS[seq(N-H+h-D,N+1-n-H+h-D,by=-1),1]
KK<-KNN.pls(X[,select.var],Y,q2[select.var],k=Kmin,C=C,F=FF,D)
p2[1:h]<-KK[(H-h+1):H]
pdirect5<-rbind(pdirect5,p2)
}
for (h in round(H/2):(H)){
p2<-numeric(H)+NA
q2<-TS[seq(N-D,N+1-n-D,by=-1),1]
KK<-KNN.pls(X[,select.var],Y[,1:h],q2[select.var],k=Kmin,C=C,F=FF,D)
p2[1:h]<-KK
pdirect5<-rbind(pdirect5,p2)
}
## combination of different predictions
p<-apply(pdirect5,2,mean,na.rm=T)
},
mimo.lin.pls={
pdirect6<-NULL
if (smooth & H >=2)
for (h in round(H/2):(H)){ ## start before the end of the series with an horizon H
p2<-numeric(H)+NA
q2<-TS[seq(N-H+h-D,N+1-n-H+h-D,by=-1),1]
KK<-lin.pls(X[,select.var],Y,q2[select.var])
p2[1:h]<-KK[(H-h+1):H]
pdirect6<-rbind(pdirect6,p2)
}
for (h in max(1,round(H-2)):(H)){
p2<-numeric(H)+NA
q2<-TS[seq(N-D,N+1-n-D,by=-1),1]
KK<-lin.pls(X[,select.var],Y[,1:h],q2[select.var])
p2[1:h]<-KK
pdirect6<-rbind(pdirect6,p2)
}
## combination of different predictions
p<-apply(pdirect6,2,mean,na.rm=T)
},
mlmimo={
pdirect3<-NULL
TS.acf<-TS[,1]
if (smooth & H >=2)
for (h in round(H/2):(H)){ ## start before the end of the series with an horizon H
p2<-numeric(H)+NA
q2<-TS[seq(N-H+h-D,N+1-n-H+h-D,by=-1),1]
KK<-pred(learner, X[,select.var],Y,q2[select.var],k=Kmin,C=C,F=FF,
TS=TS.acf,D,Reg=1)
p2[1:h]<-KK[(H-h+1):H]
pdirect3<-rbind(pdirect3,p2)
}
for (h in round(H/2):(H)){
p2<-numeric(H)+NA
q2<-TS[seq(N-D,N+1-n-D,by=-1),1]
KK<-pred(learner,X[,select.var],Y[,1:h],q2[select.var],k=Kmin,C=C,F=FF,
TS=TS.acf,D,Reg=1)
p2[1:h]<-KK
pdirect3<-rbind(pdirect3,p2)
}
## combination of different predictions
p<-apply(pdirect3,2,mean,na.rm=T)
},
iter={
piter<-numeric(H)
for (h in 1:H){
piter[h]<-KNN.multioutput(X[,select.var],array(Y[,1],c(NROW(X),1)),
q[select.var],k=Kmin,C=C,F=FF)
q<-c(piter[h],q[1:(length(q)-1)])
}
p<-piter
},
lazyiter={
piter<-numeric(H)
LPar=c(Kmin,C*Kmin)*length(select.var)
LPar[1]=min(LPar[1],NROW(X)-1)
LPar[2]=min(LPar[2],NROW(X))
if (NROW(X) <= (7*length(select.var)))
LPar=NULL
CPar=c(Kmin,C*Kmin)
CPar[1]=min(CPar[1],NROW(X)-1)
for (h in 1:H){
piter[h]<-lazy.pred(X[,select.var],array(Y[,1],c(NROW(X),1)),rbind(q[select.var]),
conPar=CPar,linPar=LPar,cmbPar=3)
q<-c(piter[h],q[1:(length(q)-1)])
}
p<-piter
},
clazyiter={
piter<-numeric(H)
LPar=NULL
CPar=c(Kmin,C*Kmin)
CPar[1]=min(CPar[1],NROW(X)-1)
for (h in 1:H){
piter[h]<-lazy.pred(X[,select.var],array(Y[,1],c(NROW(X),1)),rbind(q[select.var]),
conPar=CPar,linPar=LPar,cmbPar=3)
q<-c(piter[h],q[1:(length(q)-1)])
}
p<-piter
},
rfiter={
piter<-numeric(H)
for (h in 1:H){
piter[h]<-rf.pred(X[,select.var],array(Y[,1],c(NROW(X),1)),q[select.var],
class=FALSE,ntree=C*50)
q<-c(piter[h],q[1:(length(q)-1)])
}
p<-piter
},
liniter={
piter<-numeric(H)
for (h in 1:H){
piter[h]<-lin.pred(X[,select.var],array(Y[,1],c(NROW(X),1)),
q[select.var],class=FALSE,lambda=(1e-7)*(10*C))
q<-c(piter[h],q[1:(length(q)-1)])
}
p<-piter
},
stop("multistepAhead: Unknown method")
)
if (any(is.na(c(p))))
stop("error in multipleStepAhead")
c(p+trnd.ts) ## trend correction
}
#' MmultiplestepAhead
#' @author Gianluca Bontempi \email{Gianluca.Bontempi@@ulb.be}
#'
#' @references \url{https://tinyurl.com/sfmlh}
#' @description MmultiplestepAhead: multi-variate multi-step-ahead forecaster
#' @details Wrapper over a set of methods for multi variate multiple step ahead time series prediction
#' @title MmultiplestepAhead
#' @name MmultiplestepAhead
#' @param TS: time series [T,m] where m>1 is the number of series
#' @param n: embedding order
#' @param dfmlmethods: alternative methods from \link{MmultiplestepAhead} used by DFML to predict factors
#' @param H: horizon
#' @param unimethod: method from \link{MmultiplestepAhead} used to predict each univariate series
#' @param multi:
#' \itemize{
#' \item{UNI}: prediction based on univariate forecasting with unimethod in \link{MmultiplestepAhead}
#' \item{DFM}: prediction based on DFM
#' \item{DFML}: prediction based on DFML
#' \item{VAR}: prediction based on VAR
#' \item{VARs}: prediction based on VAR shrink from \pkg{VARshrink} package
#' \item{RNN}: prediction based on python (\pkg{reticulate}) implementation of rnn (recurrent neural network)
#' \item{LSTM}: prediction based on python (\pkg{reticulate}) implementation of lstm (long short term memory)
#' \item{TRANSF}: prediction based on python (\pkg{reticulate}) implementation of a transformer
#' \item{MIMO_rr}: prediction based on a MIMO ridge regressor (lambda estimation based on PRESS)
#' \item{MITER_rr}: prediction based on an iterated ridge regressor (lambda estimation based on PRESS and a criterion with horizon Hobj)
#' \item{MIMO_red}: prediction based on a MIMO reduced rank regressor (\pkg{rrpack})
#' \item{MIMO_cca}: prediction based on a MIMO CCA
#' \item{MIMO_pls}: prediction based on a MIMO partial least-squares (Sklearn python)
#' \item{MIMO_las}: prediction based on a MIMO lasso (Sklearn python)
#' \item{MIMO_rf}: prediction based on a MIMO Random Forest (Sklearn python)
#' \item{MIMO_ml}: prediction based on MIMO direct strategy and predictor from \link{pred}
#' \item{MIMO_fs}: prediction based on MISO direct strategy and feature selection (predictor from \link{pred})
#' }
#' @return matrix [H,m] with the H step ahead predictions for m series
#' @export
#' @details
#' The python forecasters require the installation of \pkg{reticulate} and several python packages (scikit-learn, tensorflow, keras)
#' @examples
#' ## Multi-variate Multi-step-ahead time series forecasting
#' rm(list=ls())
#' library(gbcode)
#'library(reticulate)
#'set.seed(0)
#'N=1000
#'m=3
#'n=3
#'H=10
#'TS<-array(0,c(N,m))
#'for (j in 1:m){
#' for (f in 1:5)
#' TS[,j]=TS[,j]+sin(2*pi*(1:(N))/runif(1,8,20))
#' TS[,j]=TS[,j]+rnorm(N,sd=0.3)
#'
#'}
#'TS=scale(TS)
#'N=NROW(TS)
#'
#'#P=MmultiplestepAhead(TS[1:(N-H),],n=n,H=H,multi="RNN",
#'# nepochs=100, nunits=10)
#'P=MmultiplestepAhead(TS[1:(N-H),],n=n,H=H,multi="MITER_rr",
#' nLambdas=150)
#'if (m==1)
#' P=cbind(P)
#'cat("MSE=",mean((P-TS[(N-H+1):N,])^2),"\n")
#'par(mfrow=c(1,m))
#'Nvis=round(N-4*H)
#'for (j in 1:m){
#' Yhat=numeric(N)+NA
#' Yhat[(N-H+1):N]=P[,j]
#' plot(TS[Nvis:N,j],type="l",
#' main=paste("MSE=",round(mean((TS[(N-H+1):N,j]- Yhat[(N-H+1):N])^2),2)))
#' lines(Yhat[Nvis:N],col="red",lw=3)
#'}
#'
MmultiplestepAhead<-function(TS,n=1,H=1,D=0, multi="UNI",
unimethod="stat_naive",
dfmlmodels="lindirect",
mod="rf",
pc0=2,cdfml=2,
verbose=FALSE,
debug=FALSE,...){
args<-list(...)
if (length(args)>0)
for(i in 1:length(args)) {
assign(x = names(args)[i], value = args[[i]])
}
m<-NCOL(TS)
N=NROW(TS)
if (any(is.na(TS)))
stop("Remove NA values before using this function")
if (m<=1)
stop("Only for multivariate series")
Yhat=array(NA,c(H,m))
if (multi=="UNI")
for (j in 1:m){
Yhat[,j]=multiplestepAhead(TS[,j],n,H,D=D, method=unimethod,...)
}
if (multi=="RNN")
Yhat=pyrnnpredgpt(TS,H,n,...)
if (multi=="LSTM")
Yhat=pylstmpredgpt(TS,H,n,...)
if (multi=="TRANSF")
Yhat=pytransfpredgpt(TS,H,n,...)
if (is.element(multi, c("darts_nbeats","darts_tft","darts_TCN",
"darts_transformer","darts_nhits","darts_RNN",
"darts_blockRNN","darts_lightGBM","darts_xGBM",
"darts_RandomForest" )))
Yhat=darts(TS,H,n,plearn=multi,...)
if (multi=="DFM"){
Yhat=dfml(TS,n,H,p0=pc0,dfmod=dfmlmodels[1],...)
}
if (multi=="kfm"){
Yhat=kfml(TS,n,H,p0=pc0,adaptive=FALSE,...)
}
if (multi=="kfm2"){
Yhat=kfml(TS,n,H,p0=pc0,adaptive=TRUE,...)
}
if (multi=="DFML"){
## DFML searches in the space: #Pcomponents(1:cdfml*pc0)
# #models(dfmlmodels), autoregressive order (1:cdfml*n)
Ddesign=dfmldesign(TS,cdfml*n,H,p0=cdfml*pc0,models=dfmlmodels,...)
Yhat=dfml(TS,Ddesign$n,H,p0=Ddesign$p,dfmod=Ddesign$mod,...)
}
if (multi=="VAR"){
if (n*(m^2)>N)
n=1
Yhat=VARpred2(TS,m=MTS::VAR(TS,p=n,output=FALSE),h=H,Out=FALSE)$pred
}
if (multi=="VARs"){
Yhat=VARspred(TS,n,H)
}
if (multi=="MIMO_fs")
Yhat=multifs(TS,n,H,mod=mod,debug=debug,...)
if (multi=="MIMO_rr")
Yhat=multiridge(TS,n,H,MIMO=TRUE,direct=FALSE,preq=FALSE,...)$Yhat
if (multi=="multiridge2")
Yhat=multiridge(TS,n,H,MIMO=TRUE,direct=FALSE,preq=TRUE,...)$Yhat
if (multi=="whitenridge")
Yhat=whitenridge(TS,n,H,MIMO=FALSE,direct=TRUE,...)$Yhat
if (multi=="MISO_rr")
Yhat=multiridge(TS,n,H,MIMO=FALSE,direct=TRUE,...)$Yhat
if (multi=="MIMOSO_rr")
Yhat=multiridge(TS,n,H,MIMO=TRUE,direct=TRUE,...)$Yhat
if (multi=="MITER_rr")
Yhat=multiteridge(TS,n,H,...)$Yhat
if (multi=="MITER_rr_MC")
Yhat=multiteridgeMC(TS,n,H,nLambdas=10,...)$Yhat
if (multi=="ensridge")
Yhat=ensridge(TS,n,H,...)$Yhat
if (multi=="multiridge2")
Yhat=multiridge(TS,n,H,maha=TRUE,MIMO=TRUE,...)$Yhat
if (multi=="unimultiridge"){
MR=multiridge(TS,n,H,MIMO=TRUE,...)
Yhat=MR$Yhat
for (i in 1:m){
MRi=multiridge(TS[,i],n,H,MIMO=TRUE,...)
if (MRi$MSE<MR$MSE[i])
Yhat[,i]=MRi$Yhat
}
return(Yhat)
}
if (multi=="MIMO_red")
Yhat=multirr(TS,n,H,...)
if (multi=="MIMO_pred")
Yhat=multiml(TS,n,H,learner=mod,...)
if (multi=="MIMO_las")
Yhat=multiml(TS,n,H,learner="py.lasso_regr",...)
if (multi=="MIMO_rf")
Yhat=multiml(TS,n,H,learner="py.rf_regr",...)
if (multi=="MIMO_keras")
Yhat=multiml(TS,n,H,learner="py.keras_regr",...)
if (multi=="MIMO_cca")
Yhat=multicca(TS,n,H,mod=mod,...)
if (multi=="MIMO_pls")
Yhat=multiml(TS,n,H,learner="py.pls_regr",...)
if (any(is.na(Yhat)))
stop(paste(multi," is an unknown method in MmultiplestepAhead"))
return(Yhat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.