R/predict.r

In QWDAP: Quantum Walk-Based Data Analysis and Prediction

#' @title Prediction
#' @description Based on the established model, make predict. The core algorithm of VAR prediction comes from MTS(ver. 1.1.1).
#' @usage qwdap.predict(in_model, data_range)
#'
#' @param in_model a 'QWMODEL' object, which is the model built by Stepwise Regression,
#' PCR, PLSR, PPR, VAR in this package.
#' @param data_range indicate the index range of the part data generated by quantum walks for predict.
#'
#' @return the predict data.
#' @importFrom stats predict
#' @export
#'
#' @examples
#' data(traffic.model.n1)
#' res.predict <- qwdap.predict(traffic.model.n1,c(501,720))
qwdap.predict <- function(in_model, data_range){
  if(!inherits(in_model, 'QWMODEL')){
    stop("The 'in_model' is not a 'QWMODEL' object.")
  }
  if(!is.vector(data_range)||!is.numeric(data_range)||length(data_range)<2||data_range[1]>data_range[2]){
    stop("The 'data_range' is error.")
  }
  res<-NULL
  if(in_model$method == 'VAR'){
    "VARpsi" <- function(Phi,lag=5){
      # Computes the psi-weight matrices of a VAR(p) model.
      # Phi=[phi1,phi2,phi3, ....] coefficient matrix
      # Created by Ruey S. Tsay, April 2009 & modified on April 2011.
      #
      # Compute MA representions
      k=nrow(Phi)
      m=ncol(Phi)
      p=floor(m/k)
      Si=diag(rep(1,k))
      if(p < 1) p =1
      if(lag < 1) lag=1
      #
      for (i in 1:lag){
        if (i <= p){
          idx=(i-1)*k
          tmp=Phi[,(idx+1):(idx+k)]
        }
        else{
          tmp=matrix(0,k,k)
        }
        #
        jj=i-1
        jp=min(jj,p)
        if(jp > 0){
          for(j in 1:jp){
            jdx=(j-1)*k
            idx=(i-j)*k
            w1=Phi[,(jdx+1):(jdx+k)]
            w2=Si[,(idx+1):(idx+k)]
            tmp=tmp+w1%*%w2
            ##print(tmp,digits=4)
          }
        }
        Si=cbind(Si,tmp)
      }
      
      VARpsi <- list(psi=Si)
    }
    "VARXpred" <- function(m1,newxt=NULL,hstep=1,orig=0){
      #This program predicts the VARX model.
      ## z(t) = c0 + sum_{i=1}^p phi_i * z(t-i) + \sum_{j=0}^m xt(t-j) + a(t).
      ##
      zt=as.matrix(m1$data); xt=as.matrix(m1$xt); p=m1$aror; m=m1$m
      Ph0=as.matrix(m1$Ph0); Phi=as.matrix(m1$Phi); Sig=as.matrix(m1$Sigma); beta=as.matrix(m1$beta)
      include.m=m1$include.mean
      nT=dim(zt)[1]; k=dim(zt)[2]; dx=dim(xt)[2]
      se=NULL
      if(length(Ph0) < 1)Ph0=matrix(rep(0,k),k,1)
      if(hstep < 1)hstep=1
      if(orig < 1) orig=nT
      #
      if(length(newxt) > 0){
        if(!is.matrix(newxt))newxt=as.matrix(newxt)
        ### calculate predictions
        h1=dim(newxt)[1]
        hstep=min(h1,hstep)
        nzt=as.matrix(zt[1:orig,])
        # xt=rbind(as.matrix(xt[1:orig,]),newxt)
        ### changed made on 7/30/2014
        ### changed made on 3/5/2015
        if(dx > 1){
          xt=rbind(xt[1:orig,,drop=FALSE],newxt)
        }
        else{
          xt=as.matrix(c(c(xt[1:orig,]),c(newxt)))
        }
        for (i in 1:hstep){
          tmp=Ph0
          ti=orig+i
          ## VAR part
          for (i in 1:p){
            idx=(i-1)*k
            tmp=tmp+Phi[,(idx+1):(idx+k)]%*%matrix(nzt[ti-i,],k,1)
          }
          if(m > -1){
            for (j in 0:m){
              jdx=j*dx
              tmp=tmp+beta[,(jdx+1):(jdx+dx)]%*%matrix(xt[ti-j,],dx,1)
            }
          }
          nzt=rbind(nzt,c(tmp))
        }
        ### compute standard errors of predictions
        mm=VARpsi(Phi,lag=hstep)
        Si=Sig
        se=matrix(sqrt(diag(Si)),1,k)
        if(hstep > 1){
          for (i in 2:hstep){
            idx=(i-1)*k
            wk=as.matrix(mm$psi[,(idx+1):(idx+k)])
            Si=Si+wk%*%Sig%*%t(wk)
            se1=sqrt(diag(Si))
            se=rbind(se,se1)
          }
        }
        ### Print forecasts
        cat("Point forecasts: row 'pred'. \n")
        cat("Corresponding standard errors: row 'se'. \n")
        
      }
      else{
        cat("Need new data for input variables!","\n")
      }
      #
      return <- list(pred=nzt[(orig+1):(orig+hstep),],se=se[1:hstep,],orig=orig,h=hstep)
    }
    res<-VARXpred(in_model$model, newxt = in_model$ctqw[data_range[1]:data_range[2],],
                  hstep = data_range[2]-data_range[1]+1)
  }else if(in_model$method %in% c("Stepwise Regression", "PCR", "PLSR", "PPR")){
    res<-predict(in_model$model, newdata = as.data.frame(in_model$ctqw[data_range[1]:data_range[2],]))
  }else{
    stop("This model is not supported.")
  }
  return(res)
}