########################################################################### ### An R function "test_stationarity" to test for stationarity of ### ### multivariate time series data as proposed in: ### ### ### ### Jentsch, C. and Subba Rao, S. (2015). ### ### A test for second order stationarity of a multivariate time series. ### ### Journal of Econometrics 185, No. 1, 124-161. ### ########################################################################### ### For an example see below ### ### Packages required for the function "test_stationarity" ### library(np) library(vars) test_stationarity<-function(x,Boot){ d<-dim(x)[2] n<-dim(x)[1] # maximal Number of lags used for computing the test statistic m_max<-10 # step size for the bandwidths for cross validation # for spectral density estimation h_step<-1/20 # Initialize the test statistic (with bootstrap and without bootstrap)) statistik<-rep(0,m_max+1) statistik_wb<-rep(0,m_max+1) # Defining the kernel function (Bartlett-Priestley) Kern<-function(h,x){ Kern<-matrix(0,length(h),length(x)) for(j in 1:length(h)){ for(k in 1:length(x)){ if(abs(x[k])<=pi*h[j]){Kern[j,k]<-(2*pi)*1/h[j]*(3/(4*pi))*(1-(x[k]/(pi*h[j]))^2)} else{Kern[j,k]<-0} } } Kern } # Cholesky decomposition for a complex matrix cchol<-function(x){ cchol<-matrix(rep(0,(dim(x)[1]*dim(x)[1])),dim(x)[1],dim(x)[2]) for(k in 1:(dim(x)[1])){ for(l in 1:k){ hilf<-0 if ((k-1)>0 && k==l) hilf<-t(cchol[k,(1:(k-1))])%*%Conj(cchol[k,(1:(k-1))]) if (k==l) cchol[k,l]<-sqrt(x[k,l]-hilf) if ((l-1)>0 && k>l) hilf<-t(cchol[k,(1:(l-1))])%*%Conj(cchol[l,(1:(l-1))]) if (k>l) cchol[k,l]<-(1/cchol[l,l])*(x[k,l]-hilf) } } cchol } ################################################ ### spectral density estimation ### ################################################ # Computing the matrix of mDFT of data x at frequencies # \omega_0,...,\omega_{n-1} Jn<-mvfft(x)/sqrt(2*pi*n) #Jn # Arranging a bigger matrix of mDFT of data x at frequencies # \omega_{-n},...,\omega_{2n-1} Jnbig<-matrix(rep(0,3*n*d),(3*n),d) Jnbig<-rbind(Jn[1,],Conj(Jn[n:1,]),Jn[2:n,],Jn[1,],Conj(Jn[n:2,])) #Jnbig # Arranging fourier frequencies -pi to 2*pi*(2*n-1)/n freq<-rep(0,3*n) freq<-(2*pi/n)*(-n):(2*n-1) # Arranging a large matrix containing the periodograms Inbig<-matrix(rep(0,3*n*d^2),d,3*d*n) for(j in 1:(3*n)){ Inbig[,((j-1)*d+1):(j*d)]<-Jnbig[j,]%*%t(Conj(Jnbig[j,])) } ########################################################################## ### Computing the bandwidth h_CV using CV criteria ### ### on a grid of [0,1] with step size h_step ### ########################################################################## # g starts at ((2/(h_step*n))+1) not to run into singularities CV<-rep(0,2/h_step) CV_new<-rep(0,2/h_step) for(g in 2:(2/h_step)){ h<-g*h_step omega<-0 kernvec<-rep(0,n) fhut_mult_cv<-t(matrix(rep(0,n*d*d/2),(n*d)/2,d)) #fhut_mult<-t(matrix(rep(0,n*d*d/2),(n*d)/2,d)) #rep(0,n/2) for(k in 1:floor(n/2)){ omega<-(2*pi/n)*k kernvec<-Kern(h,omega-freq) kernvec[k+1]<-0 kernvec[n+1-k]<-0 kernvec[n+1+k]<-0 kernvec[2*n+1-k]<-0 kernvec[2*n+1+k]<-0 kernvec[3*n+1-k]<-0 kernvec_matrix<-kronecker(kernvec,diag(1,d)) # The conjugate due to transposing Inbig (nein) #fhut_mult_cv[,((k-1)*d+1):(k*d)]<-(1/n)*(kernvec_matrix%*%t(Inbig)) fhut_mult_cv[,((k-1)*d+1):(k*d)]<-(1/n)*(Inbig%*%t(kernvec_matrix)) } #fhut_mult_cv[,((k-1)*d+1):(k*d)] blub<-0 for(j in 1:floor(n/2)){ #Achtung hier die Indizes so, dass fhut_mult_cv und Inbig passen blub<-log(prod(eigen(fhut_mult_cv[,((j-1)*d+1):(j*d)])$values))+sum(diag(Inbig[,((n+j+1-1)*d+1):((n+j+1)*d)]%*%solve(fhut_mult_cv[,((j-1)*d+1):(j*d)])))+blub } CV_new[g]<-blub } #CV_new CV_new<-Re(CV_new) # Choosing the minimizer ((+1 because the loop over g starts in second entry) h_CV_new<-(which.min(CV_new[2:(2/h_step)])+1)*h_step #h.vector[t]<-h_CV_new h<-h_CV_new #h ####################################################################################################### # spectral density matrix estimation at # fourier frequencies omega_1,...,omega_n (more efficient implementation!) omega<-0 kernvec<-rep(0,n) fhut_mult<-t(matrix(rep(0,n*d*d),(n*d),d)) for(k in 1:n){ omega<-(2*pi/n)*k kernvec_matrix<-kronecker(Kern(h,omega-freq),diag(1,d)) fhut_mult[,((k-1)*d+1):(k*d)]<-(1/n)*(Inbig%*%t(kernvec_matrix)) ###fhut_mult[,((k-1)*d+1):(k*d)]<-(1/n)*(kernvec_matrix%*%t(Inbig)) } #fhut_mult # Rearranging and transposing Jn to get Jn_help at frequencies # omega_1,...,omega_n Jn_help<-t(rbind(Jn[2:n,],Jn[1,])) # Computing the standardized DFT Jn_stand Jn_stand<-matrix(rep(0,n*d),d,n) for(k in 1:n){ Jn_stand[,k]<-solve(cchol(fhut_mult[,((k-1)*d+1):(k*d)]))%*%Jn_help[,k] } # Computing covariances of the DFTs from lag 0 to lag n-1 Jn_stand_long<-cbind(Jn_stand,Jn_stand) DFTcov<-matrix(rep(0,(m_max+1)*d*d),d,(m_max+1)*d) for(j in 1:(m_max+1)){ for(k in 1:d){ for(i in 1:d){ DFTcov[k,((j-1)*d+i)]<-(1/n)*Jn_stand_long[k,1:n]%*%Conj(Jn_stand_long[i,((1+j-1):(n+j-1))]) } } } # Creating a matrix that contains the vectorized DFT covariances DFTcovvec<-matrix(rep(0,n*d*d),d*d,n) for(k in 1:(m_max+1)){ blub<-DFTcov[,((k-1)*d+1):(k*d)] dim(blub)<-NULL DFTcovvec[,k]<-blub } ######################################################################### ### Stationary Bootstrap begins here ### ######################################################################### ######################################################################## ### Rule-of-thumb for the tuning parameter l (Politis & White (2003) ### ######################################################################## # Block length selection s<-b.star(x,round=FALSE)[,1] # taking the average of the component-wise chosen block length l<-mean(s) #s #l p<-1/l #p ######################################################################## ### Stationary bootstrap loop begins here ### ######################################################################## DFTcovstern_all<-matrix(rep(0,(m_max+1)*d*d*Boot),d*d*(m_max+1),Boot) for(b in 1:Boot){ # Generating the bootstrap series xstern<-matrix(rep(0,10*n*d),10*n,d) xlong<-rbind(x,x,x,x,x,x,x,x,x,x) length<-0 while(lengthsqrt(2.4*log(n))){ vector_penalty<-seq(1,m_max,by=1)*2 } if(max(sqrt(n)*abs(DFTcovvech_stand_acc_wb[2:(m_max+1)]))>sqrt(2.4*log(n))){ vector_penalty_wb<-seq(1,m_max,by=1)*2 } # Setting the chosen number of lags used for the test statistic (with bootstrap) m_opt<-1 maxima<-which((vector_portmanteau[2:m_max]-vector_penalty[2:m_max])-(vector_portmanteau[1:m_max-1]-vector_penalty[1:m_max-1])<=0) m_opt<-maxima[1] if(is.na(maxima[1])==TRUE) m_opt<-m_max # Setting the chosen number of lags used for the test statistic (without bootstrap) m_opt_wb<-1 maxima_wb<-which((vector_portmanteau_wb[2:m_max]-vector_penalty_wb[2:m_max])-(vector_portmanteau_wb[1:m_max-1]-vector_penalty_wb[1:m_max-1])<=0) m_opt_wb<-maxima_wb[1] if(is.na(maxima_wb[1])==TRUE) m_opt_wb<-m_max ############################################## ### Here ends the selection procedure of m ### ############################################## # Computing the test statistics based on DFT covariance lags up to 1,2,3,...,m_max for(m in 1:m_max){ for(k in 1:m){ statistik[m]<-sum((DFTcovvech_stand_re[,k+1])^2)+sum((DFTcovvech_stand_im[,k+1])^2)+statistik[m] statistik_wb[m]<-sum((DFTcovvech_stand_re_wb[,k+1])^2)+sum((DFTcovvech_stand_im_wb[,k+1])^2)+statistik_wb[m] } } # Computing the test statistics with automatically chosen number of DFT covariance lags for(k in 1:m_opt){ statistik[m_max+1]<-sum((DFTcovvech_stand_re[,k+1])^2)+sum((DFTcovvech_stand_im[,k+1])^2)+statistik[m_max+1] statistik_wb[m_max+1]<-sum((DFTcovvech_stand_re_wb[,k+1])^2)+sum((DFTcovvech_stand_im_wb[,k+1])^2)+statistik_wb[m_max+1] } cbind(statistik,qchisq(1-0.05,c(1:m_max,m_opt)*d*(d+1)),statistik_wb,qchisq(1-0.05,c(1:m_max,m_opt_wb)*d*(d+1))) } ################################################ ### End of Definition of "test_stationarity" ### ################################################ ################## ### An Example ### ################## library(MASS) ###Parameters to create a bivariate AR(1) sample ### # Dimension of the process d<-3 # sample size of multivariate (d-dimensional) data n<-100 # Parameters for VAR(1) model if(d==2){ A<-matrix(c(0.5,0,-0.3,0.7),2,2) Sigma<-matrix(c(1,0,0,1),2,2) } if(d==3){ A<-matrix(c(0.5,0,0,0,-0.3,0.7,0,0.2,0.3),3,3) Sigma<-matrix(c(1,0,0,0,1,0,0,0,1),3,3) } # Generate the time series sample x<-matrix(0,d,n+100) e<-mvrnorm(n+100,rep(0,d),Sigma,d,d) e<-t(e) for(k in 2:(n+100)){ x[,k]<-A%*%x[,(k-1)]+e[,k] } x<-x[,101:(100+n)] x<-t(x) ### the time series data x to be testet shall be a n \times d matrix ### ### Boot is the number of bootstrap replications ### Boot<-250 test_stationarity(x,Boot) # gives the test statistic based on DFT covariance up to lag 1,2,3,...,m_max # and with automatically chosen optimal m in the last row together with # corresponding critical values from chi^2 distribution at 95% level # with stationary bootstrap variance adjustment (first two columns) # and without bootstrap (Gaussian case) (last two columns)