Nothing
#' @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(!inherits(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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.