R/var.r

In QWDAP: Quantum Walk-Based Data Analysis and Prediction

#' @title Vector Autoregressive Model
#' @description Vector autoregressive model. This is a regression method used to
#' establish the temporal relationship between the original time series and
#' the modes generated by quantum walks. The core algorithm comes from MTS(ver. 1.1.1).
#' @usage qwdap.var(in_data, data_range, plotting)
#' @param in_data a 'QWMS' object, which includes the target series and the
#' selected modes which can be obtained from modes selection.
#' @param data_range the range of the train samples.
#' @param plotting whether to plot.
#'
#' @return a 'QWMODEL' object which includes the information of regression analysis.
#' @importFrom graphics lines legend
#' @export qwdap.var
#'
#' @examples
#' data("traffic.n1")
#' res.var <- qwdap.var(traffic.n1,c(1,500))
#' 
qwdap.var<-function(in_data, data_range, plotting = FALSE){
  if(class(in_data)!='QWMS'){
    stop("The 'in_data' is not a 'QWMS' object.")
  }
  if(!is.vector(data_range)||!is.numeric(data_range)||length(data_range)<2){
    stop("The parameter 'data_range' is error.")
  }
  # order
  "VARXorder" <- function(x,exog,maxp=13,maxm=3,output=F){
    # Compute the AIC, BIC, HQ values and M-stat
    ##### This is a modified version of the old program in "VARorderE",
    ##### which uses the same number of data points.
    #####
    x1=as.matrix(x)
    exog=as.matrix(exog)
    nT=dim(x1)[1]
    k=dim(x1)[2]
    ksq=k*k
    if(maxp < 1)maxp=1
    nT1=dim(exog)[1]; m=dim(exog)[2]
    #
    if(nT1 > nT){
      cat("Adjustment made for different nobs:",c(nT,nT1), "\n")
    }
    if(nT > nT1){
      cat("Adjustment made for different nobs:",c(nT,nT1),"\n")
      nT=nT1
    }
    ###
    aic=matrix(0,maxp+1,maxm+1)
    bic=aic; hq=aic
    for (mm in 0:maxm){
      ### start with VAR(0) model, which uses just the sample means.
      isto=mm+1
      y=x1[isto:nT,]
      xm=rep(1,rep(nT-mm))
      for (i in 0:mm){
        xm=cbind(xm,exog[(isto-i):(nT-i),])
      }
      xm=as.matrix(xm)
      xpx=crossprod(xm,xm)
      xpxi=solve(xpx)
      xpy=t(xm)%*%y
      beta=xpxi%*%xpy
      y=y-xm%*%beta
      #
      s=t(y)%*%y/(nT-mm)
      enob=nT-mm
      c1=log(det(s))
      aic[1,mm+1]=c1+2*k*mm/enob
      bic[1,mm+1]=c1+log(enob)*k*m/enob
      hq[1,mm+1]=c1+2*log(log(enob))*k*m/enob
      #
      for (p in 1:maxp){
        ist=max(mm+1,p+1)
        enob=nT-ist+1
        y=as.matrix(x1[ist:nT,])
        xmtx=rep(1,enob)
        for (i in 0:mm){
          xmtx=cbind(xmtx,exog[(ist-i):(nT-i),])
        }
        for (j in 1:p){
          xmtx=cbind(xmtx,x1[(ist-j):(nT-j),])
        }
        xm1=as.matrix(xmtx)
        xpx = t(xm1)%*%xm1
        xpxinv=solve(xpx)
        xpy=t(xm1)%*%y
        beta=xpxinv%*%xpy
        resi=y-xm1%*%beta
        sse=(t(resi)%*%resi)/enob
        #print(paste("For p = ",p,"residual variance is", sse))
        d1=log(det(sse))
        npar=p*ksq+k*(mm+1)
        aic[p+1,mm+1]=d1+(2*npar)/enob
        bic[p+1,mm+1]=d1+(log(enob)*npar)/enob
        hq[p+1,mm+1]=d1+(2*log(log(enob))*npar)/enob
      }
      #end of for (mm in 0:maxm)
    }
    
    ind.min <- function(x){
      if(!is.matrix(x))x=as.matrix(x)
      r=dim(x)[1]
      c=dim(x)[2]
      xm=min(x)
      COntin=TRUE
      while(COntin){
        for(j in 1:c){
          idx=c(1:r)[x[,j]==xm]
          jj=j
          if(length(idx) > 0){
            ii=c(1:r)[x[,jj]==xm][1]
            kdx=c(ii,jj)
            ##print(kdx)
            COntin=FALSE
          }
        }
      }
      kdx
    }
    ## selection
    aicor=ind.min(aic)-1; bicor=ind.min(bic)-1;hqor=ind.min(hq)-1
    if(output){
      cat("selected order(p,s): aic = ",aicor,"\n")
      cat("selected order(p,s): bic = ",bicor,"\n")
      cat("selected order(p,s): hq = ",hqor,"\n")
    }
    
    VARXorder<-list(aic=aic,aicor=aicor,bic=bic,bicor=bicor,hq=hq,hqor=hqor)
  }
  # var
  "VARX" <- function(zt,p,xt=NULL,m=0,include.mean=T,fixed=NULL,output=T){
    #This command fits the 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(zt)
    if(length(xt) < 1){
      m = -1; kx=0}
    else{
      xt=as.matrix(xt); kx=dim(xt)[2]
    }
    if(p < 0)p=0
    ist=max(p,m)+1
    nT=dim(zt)[1]
    k=dim(zt)[2]
    yt=zt[ist:nT,]
    xmtx=NULL
    if(include.mean)xmtx=rep(1,(nT-ist+1))
    #
    if(p > 0){
      for (i in 1:p){
        xmtx=cbind(xmtx,zt[(ist-i):(nT-i),])
      }
    }
    #
    if( m > -1){
      for (j in 0:m){
        xmtx=cbind(xmtx,xt[(ist-j):(nT-j),])
      }
    }
    #
    p1=dim(xmtx)[2]
    nobe=dim(xmtx)[1]
    ##cat("dim of xmtx",c(nobe,p1),"\n")
    #
    if(length(fixed) < 1){
      ## no constriants
      xpx=t(xmtx)%*%xmtx
      xpy=t(xmtx)%*%yt
      xpxi=solve(xpx)
      beta=xpxi%*%xpy
      resi=as.matrix(yt-xmtx%*%beta)
      sig=crossprod(resi,resi)/nobe
      co=kronecker(sig,xpxi)
      se=sqrt(diag(co))
      se.beta=matrix(se,nrow(beta),k)
      npar=nrow(beta)*k
      d1=log(det(sig))
      aic=d1+2*npar/nobe
      bic=d1+(npar*log(nobe))/nobe
    }
    else{
      # with zero-parameter constriants
      beta=matrix(0,p1,k)
      se.beta=matrix(1,p1,k)
      resi=yt
      npar=0
      for (i in 1:k){
        idx=c(1:p1)[fixed[,i] > 0]
        npar=npar+length(idx)
        if(length(idx) > 0){
          xm=as.matrix(xmtx[,idx])
          y1=matrix(yt[,i],nobe,1)
          xpx=t(xm)%*%xm
          xpy=t(xm)%*%y1
          xpxi=solve(xpx)
          beta1=xpxi%*%xpy
          res = y1 - xm%*%beta1
          sig1=sum(res^2)/nobe
          se=sqrt(diag(xpxi)*sig1)
          beta[idx,i]=beta1
          se.beta[idx,i]=se
          resi[,i]=res
        }
        # end of for (i in 1:k)
      }
      #
      sig=crossprod(resi,resi)/nobe
      d1=log(det(sig))
      aic=d1+2*npar/nobe
      bic=d1+log(nobe)*npar/nobe
      # end of else
    }
    ### print
    Ph0=NULL
    icnt=0
    if(include.mean){
      Ph0=beta[1,]; icnt=icnt+1
      if(output){
        cat("constant term: ","\n")
        cat("est: ",round(Ph0,4),"\n")
        cat(" se: ",round(se.beta[1,],4),"\n\n")
      }
    }
    Phi=NULL
    if(p > 0){
      Phi=t(beta[(icnt+1):(icnt+k*p),])
      sePhi=t(se.beta[(icnt+1):(icnt+k*p),])
      for (j in 1:p){
        jcnt=(j-1)*k
        if(output){
          cat("AR(",j,") matrix")
          print(round(Phi[,(jcnt+1):(jcnt+k)],3))
          cat("standard errors","\n")
          print(round(sePhi[,(jcnt+1):(jcnt+k)],3))
          cat("\n")
        }
      }
      icnt=icnt+k*p
      ## end of if(p > 0)
    }
    if(m > -1){
      if(output){
        cat("Coefficients of exogenous","\n")
      }
      Beta=t(beta[(icnt+1):(icnt+(m+1)*kx),])
      seBeta=t(se.beta[(icnt+1):(icnt+(m+1)*kx),])
      if(kx == 1){
        Beta=t(Beta)
        seBeta=t(seBeta)
      }
      for (i in 0:m){
        jdx=i*kx
        if(output){
          cat("lag-",i," coefficient matrix","\n")
          print(round(Beta[,(jdx+1):(jdx+kx)],3))
          cat("standard errors","\n")
          print(round(seBeta[,(jdx+1):(jdx+kx)],3))
        }
      }
      ## end of if(m > -1)
    }
    ##
    if(output){
      cat("Residual Covariance Matrix","\n")
      print(round(sig,5))
      cat("===========","\n")
      cat("Information criteria: ","\n")
      cat("AIC: ",aic,"\n")
      cat("BIC: ",bic,"\n")
    }
    VARX <- list(data=zt,xt=xt,aror=p,m=m,Ph0=Ph0,Phi=Phi,beta=Beta,residuals=resi,Sigma=sig,
                 coef=beta,se.coef=se.beta,include.mean=include.mean)
  }
  # pre combine
  co_data <- cbind(in_data$real, in_data$ctqw)
  co_data <- subset(co_data, select = c(colnames(in_data$real), in_data$variate))
  indexes<-VARXorder(co_data[data_range[1]:data_range[2],1],
                     co_data[data_range[1]:data_range[2],-1],10)
  if(indexes$aicor[1]==0){
    indexes$aicor[1]=indexes$bicor[1]
  }
  ## 使用VARX
  res <- VARX(co_data[data_range[1]:data_range[2],1],
              indexes$aicor[1],co_data[data_range[1]:data_range[2],-1],
              indexes$aicor[2],output = F)
  if(plotting){
    plot(x=c(data_range[1]:data_range[2]),
         y=co_data[data_range[1]:data_range[2],1],
         type="l",col=1,xlab="index",ylab="value",lwd=1)
    cir<-c((data_range[2]-length(res$residuals[,1])):(data_range[2]-1))
    lines(x=cir,y=res$residuals[,1]+co_data[cir,1],type="l",col=2,lwd=1)
    legend("topleft", c("Actual series","Fitted series"), col = c(1,2),
           lwd = c(1), bg = "grey95", box.col = NA,
           cex = 0.8, inset = c(0.02, 0.03), ncol = 1)
  }
  res<-list(real = in_data$real, ctqw = co_data[,-1], index = in_data$index,
            method = "VAR",model=res)
  res<-structure(res,class="QWMODEL")
  return(res)
}