################################################################# ## X is the data matrix ## ## ## ## Boot is the number of bootstrap samples to be generated ## ## ## ## l is the banding parameter; default is l=1 ## ## ## ## eps and beta are the parameters for making ## ## Gamma_kappa_matrix positive definite if necessary; ## ## default value are eps=1 and beta=1 ## ## ## ## if l_automatic==1 the banding parameter is chosen ## ## data-adaptively; default value is l_automatic==0 ## ## ## ## if l_automatic_local==0 the banding parameter is chosen ## ## globally and if l_automatic_local==1 locally for each ## ## covariance component; default value is l_automatic_local==0 ## ################################################################# MLPB<-function(X,Boot,l=1,eps=1,beta=1,l_automatic=0,l_automatic_local=0){ # Dimension of observed values n1<-dim(X)[1] # sample sizes of (n1-dimensional) data observed n2<-dim(X)[2] N<-n1*n2 # l_matrix<-matrix(rep(l,n1^2),n1,n1) # Defining function kappa (trapezoid [cf. McMurry & Politis (2010)]) kappa<-function(x,l){ kappa<-matrix(0,dim(l)[1],dim(l)[2]) for(j in 1:dim(l)[1]){ for(k in 1:dim(l)[2]){ if(l[j,k]>0){ kappa[j,k]<-2-abs(x/l[j,k]) if(abs(x/l[j,k])<=1) kappa[j,k]<-1 if(abs(x/l[j,k])>2) kappa[j,k]<-0 } } } kappa } # Rearranging the data columnwise in one large n1*n2-dimensional vector X_vec<-X dim(X_vec)<-NULL ######################################################################## ### Here starts the computation of covariances ### ######################################################################## # Estimating all multivariate covariances and putting them in # one large n1*(2n2-1) matrix C_matrix<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1) C_matrix_kappa_trans<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1) R_matrix_kappa_trans<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1) # Sample mean if(n1==1) X_mean<-mean(X) if(n1!=1){ X_mean<-rep(0,n1) for(k in 1:n1){ X_mean[k]<-mean(X[k,]) } } # Computing C_matrix for(k in (-(n2-1)):(n2-1)){ blub<-0 if(n1==1) C_matrix[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-(1/n2)*t(X[(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%(X[(max(1,1-k)):(min(n2,n2-k))]-X_mean) if(n1!=1) C_matrix[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-(1/n2)*(X[,(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%t(X[,(max(1,1-k)):(min(n2,n2-k))]-X_mean) } if(l_automatic==1){ ################################## ### Computing Autocorrelations ### ################################## R_matrix<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1) denominator_matrix<-matrix(rep(0,n1*n1),n1,n1) diag(denominator_matrix)<-sqrt(diag(C_matrix[,(n1*(n2-1)+1):(n1*n2)])) denominator_matrix<-solve(denominator_matrix) for(k in 1:(2*n2-1)){ R_matrix[,((k-1)*n1+1):(k*n1)]<-denominator_matrix%*%C_matrix[,((k-1)*n1+1):(k*n1)]%*%denominator_matrix } if(l_automatic_local==0){ M_zero<-2 h<-1 marker<-0 while(marker==0){ #Kn_vec<-seq(1,max(5,log10(n2))) #if(max((abs(R_matrix[,((2*n2-1)+h*n1):((2*n2-1+n1-1)+n1*(max(5,log10(n2))+h-1))])))0) C_matrix_kappa_trans[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-t((kappa(k,l_matrix)/n2)*(X[,(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%t(X[,(max(1,1-k)):(min(n2,n2-k))]-X_mean)) } # Computing the diagonal matrix for standardization to get autocorrelations # (as above) denominator_matrix<-matrix(rep(0,n1*n1),n1,n1) diag(denominator_matrix)<-sqrt(diag(C_matrix[,(n1*(n2-1)+1):(n1*n2)])) denominator_matrix<-solve(denominator_matrix) # Computing R_matrix_kappa_trans (for equivariance to construct positive # definite estimator) for(k in (-(n2-1)):(n2-1)){ blub<-0 if(n1==1) R_matrix_kappa_trans[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-(1/var(X))*(kappa(k,l_matrix)/n2)*t(X[(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%(X[(max(1,1-k)):(min(n2,n2-k))]-X_mean) if(n1!=1 && k<0) R_matrix_kappa_trans[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-denominator_matrix%*%t(t((kappa(k,l_matrix)/n2))*(X[,(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%t(X[,(max(1,1-k)):(min(n2,n2-k))]-X_mean))%*%denominator_matrix if(n1!=1 && k==0) R_matrix_kappa_trans[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-denominator_matrix%*%t((1/n2)*(X[,(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%t(X[,(max(1,1-k)):(min(n2,n2-k))]-X_mean))%*%denominator_matrix if(n1!=1 && k>0) R_matrix_kappa_trans[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-denominator_matrix%*%t((kappa(k,l_matrix)/n2)*(X[,(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%t(X[,(max(1,1-k)):(min(n2,n2-k))]-X_mean))%*%denominator_matrix } # Arranging the autocovariance matrices in R_matrix_kappa_trans in one large # matrix corresponding to X_vec R_kappa_matrix<-matrix(rep(0,N*N),N,N) for(i in 1:n2){ R_kappa_matrix[(((i-1)*n1+1):((i-1)*n1+n1)),]<-R_matrix_kappa_trans[,(n1*(n2-i)+1):(n1*(2*n2-i))] } R_kappa_matrix_hilf<-R_kappa_matrix # Computing the eigenvalues and eigenvectors of Gamma_kappa_matrix spectraldecomposition<-eigen(R_kappa_matrix) # orthogonal matrix S<-spectraldecomposition$vectors # diagonal matrix D<-diag(spectraldecomposition$values) # Manipulating the diagonal matrix D to ensure positive definiteness if(min(diag(D))<=eps*N^(-beta)){ D_eps_counter<-D_eps_counter+1 D_eps<-D D_non_definite<-D for(k in 1:N){ D_eps[k,k]<-max(D[k,k],eps*N^(-beta)) } # Computing the banded positive definite covariance matrix estimator R_kappa_matrix<-S%*%D_eps%*%t(S) } # Get back to autocovariances! Gamma_kappa_matrix<-matrix(rep(0,N*N),N,N) nominator_matrix<-solve(denominator_matrix) for(i in 1:n2){ for(j in 1:n2){ Gamma_kappa_matrix[(((i-1)*n1+1):((i-1)*n1+n1)),(((j-1)*n1+1):((j-1)*n1+n1))]<-nominator_matrix%*%R_kappa_matrix[(((i-1)*n1+1):((i-1)*n1+n1)),(((j-1)*n1+1):((j-1)*n1+n1))]%*%nominator_matrix } } # Computing the (lower-left) Cholesky decomposition L<-t(chol(Gamma_kappa_matrix)) L_inv<-solve(L) # Multilpication of L and centered X_vec that is X_vec_cent X_cent<-matrix(rep(0,n1*n2),n1,n2) for(k in 1:n1){ X_cent[k,]<-X[k,]-X_mean[k] } # Rearranging the data columnwise in one large n1*n2-dimensional vector X_vec_cent<-X_cent dim(X_vec_cent)<-NULL # "Whitening" the centered data W_vec<-L_inv%*%X_vec_cent # Computing the centered and standardized residuals Z_vec W_vec_cent<-W_vec-mean(W_vec) W_vec_var<-var(W_vec) dim(W_vec_var)<-NULL Z_vec<-(1/sqrt(W_vec_var))*W_vec_cent W_matrix<-matrix(W_vec,n1,n2) W_matrix_cent<-matrix(rep(0,n1*n2),n1,n2) for(k in 1:n1){ W_matrix_cent[k,]<-W_matrix[k,]-mean(W_matrix[k,]) } W_matrix_var<-(1/n2)*W_matrix_cent%*%t(W_matrix_cent) W_matrix_var_inv<-solve(t(chol(W_matrix_var))) W_matrix_stud<-matrix(rep(0,n1*n2),n1,n2) for(k in 1:n2){ W_matrix_stud[,k]<-W_matrix_var_inv%*%W_matrix_cent[,k] } Z_vec<-W_matrix_stud dim(Z_vec)<-NULL ########################################################################## ### Here starts the LPB bootstrap loop ### ########################################################################## X_stern_all<-matrix(0,n1*Boot,n2) for(b in 1:Boot){ # Resampling sub-vectors from Z_vec Z_vec_stern<-rep(0,n1*n2) sv_draw<-sample(1:n2,n2,replace=TRUE) for(k in 1:n2){ Z_vec_stern[((k-1)*n1+1):(k*n1)]<-Z_vec[((sv_draw[k]-1)*n1+1):(sv_draw[k]*n1)] } # Computing the bootstrap (vector-) sample X_vec_stern X_vec_stern<-L%*%Z_vec_stern X_vec_stern # Rearranging the data vector in a data matrix X_stern<-matrix(rep(0,n1*n2),n1,n2) for(j in 1:n2){ X_stern[,j]<-X_vec_stern[((j-1)*n1+1):(j*n1)] } X_stern_all[((b-1)*n1+1):((b*n1)),]<-X_stern } ########################################################################## ### Here ends the LPB bootstrap loop ### ########################################################################## X_stern_all }