R/CP_functions_unified.R In HDTSA: High Dimensional Time Series Analysis Tools

Documented in CP_MTSDGP.CP

#' @title Estimation of matrix CP-factor model
#' @description \code{CP_MTS()} deals with CP-decomposition for high-dimensional
#'  matrix time series proposed in Chang et al. (2023):\deqn{{\bf{Y}}_t = {\bf A \bf X}_t{\bf B}^{'} +
#' {\boldsymbol{\epsilon}}_t, } where \eqn{{\bf X}_t = diag(x_{t,1},\ldots,x_{t,d})} is an \eqn{d \times d}
#' latent process, \eqn{{\bf A}} and \eqn{{\bf B}} are , respectively, \eqn{p
#' \times d} and \eqn{q \times d} unknown constant matrix, and \eqn{ {\boldsymbol{\epsilon}}_t }
#'  is a \eqn{p \times q} matrix white noise process. This function aims to estimate the rank
#'  \eqn{d} and the coefficient matrices \eqn{{\bf A}} and \eqn{{\bf B}}.
#'
#' @param Y A \eqn{n \times p \times q} data array, where \eqn{n} is the sample size and \eqn{(p,q)}
#' is the dimension of \eqn{{\bf Y}_t}.
#' @param xi A \eqn{n \times 1} vector. If \code{NULL} (the default), then a PCA-based \eqn{\xi_{t}}
#' is used [See Section 5.1 in Chang et al. (2023)] to calculate the sample auto-covariance matrix
#' \eqn{\widehat{\bf \Sigma}_{\bf Y, \xi}(k)}.
#' @param Rank A list of the rank \eqn{d},\eqn{d_1} and \eqn{d_2}. Default to \code{NULL}.
#' @param lag.k Integer. Time lag \eqn{K} is only used in \code{CP.Refined} and \code{CP.Unified} to
#' calculate the nonnegative definte matrices \eqn{\widehat{\mathbf{M}}_1} and
#' \eqn{\widehat{\mathbf{M}}_2}: \deqn{\widehat{\mathbf{M}}_1\ =\
#'   \sum_{k=1}^{K}\widehat{\mathbf{\Sigma}}_{\bf Y, \xi}(k)\widehat{\mathbf{\Sigma}}_{\bf Y, \xi}(k)',
#'   }, \deqn{\widehat{\mathbf{M}}_2\ =\
#'   \sum_{k=1}^{K}\widehat{\mathbf{\Sigma}}_{\bf Y, \xi}(k)'\widehat{\mathbf{\Sigma}}_{\bf Y, \xi}(k),
#'   }
#'   where \eqn{\widehat{\mathbf{\Sigma}}_{\bf Y, \xi}(k)} is the sample auto-covariance of
#'   \eqn{ {\bf Y}_t} and \eqn{\xi_t} at lag \eqn{k}.
#' @param lag.ktilde Integer. Time lag \eqn{\tilde K} is only used in \code{CP.Unified} to calulate the
#' nonnegative definte matrix \eqn{\widehat{\mathbf{M}}}: \deqn{\widehat{\mathbf{M}} \ =\
#'   \sum_{k=1}^{\tilde K}\widehat{\mathbf{\Sigma}}_{\tilde{\bf Z}}(k)\widehat{\mathbf{\Sigma}}_{\tilde{\bf Z}}(k)'.
#'   }
#' @param method Method to use: \code{CP.Direct} and \code{CP.Refined}, Chang et al.(2023)'s direct and refined estimators;
#'  \code{CP.Unified}, Chang et al.(2024+)'s unified estimation procedure.
#'
#' @return An object of class "mtscp" is a list containing the following
#'   components:
#'   \item{f}{The estimated latent process \eqn{(\hat{x}_{1,t},\ldots,\hat{x}_{d,t})}.}
#'   \item{Rank}{The estimated rank \eqn{(\hat{d}_1,\hat{d}_2,\hat{d})} of the matrix CP-factor model.}
#'
#'
#' @references
#'   Chang, J., He, J., Yang, L. and Yao, Q.(2023). \emph{Modelling matrix time series via a tensor CP-decomposition}.
#'   Journal of the Royal Statistical Society Series B: Statistical Methodology, Vol. 85(1), pp.127--148.
#'
#'   Chang, J., Du, Y., Huang, G. and Yao, Q.(2024+). \emph{On the Identification and Unified Estimation
#'   Procedure for the Matrix CP-factor Model}, Working paper.
#'
#' @examples
#' p = 10
#' q = 10
#' n = 400
#' d = d1 = d2 = 3
#' data <- DGP.CP(n,p,q,d1,d2,d)
#' Y = data$Y #' res1 <- CP_MTS(Y,method = "CP.Direct") #' res2 <- CP_MTS(Y,method = "CP.Refined") #' res3 <- CP_MTS(Y,method = "CP.Unified") #' @export #' @useDynLib HDTSA #' @importFrom stats arima.sim rnorm runif CP_MTS = function(Y,xi = NULL, Rank = NULL, lag.k = 15, lag.ktilde = 10, method = c("CP.Direct","CP.Refined","CP.Unified")){ n = dim(Y)[1]; p = dim(Y)[2]; q = dim(Y)[3]; if(is.null(xi)){ xi = est.xi(Y) } method <- match.arg(method) if(method == "CP.Direct"){ S_yxi_1 = Autocov_xi_Y(Y,xi,lag.k = 1) S_yxi_2 = Autocov_xi_Y(Y,xi,lag.k = 2) if(p > q){ ##(1) estimation of d K1 = t(S_yxi_1)%*%S_yxi_1 eg1 = eigen(K1) w = eg1$values
ww = w[-1]/w[-length(w)]
d = which(ww==min(ww[1:floor(0.75*q)]))

if (d > 1){
K1 = eg1$vectors[,1:d]%*%diag(eg1$values[1:d])%*%t(eg1$vectors[,1:d]); }else{ K1 = eg1$vectors[,1]%*%diag(eg1$values[1],1)%*%t(eg1$vectors[,1]);
}
K2 = t(S_yxi_1)%*%S_yxi_2;

##(2) estimation of A and B
Geg = geigen::geigen(K2,K1);
evalues = Geg$values[which(Mod(Geg$values)<=10^5&Geg$values!=0)] Bl = Geg$vectors[,which(Geg$values %in% evalues)] A = apply(S_yxi_1%*%Bl,2,l2s) Al = t(MASS::ginv(A)) B = apply(t(S_yxi_1)%*%Al,2,l2s) }else{ ##(1) estimation of d K1 = S_yxi_1%*%t(S_yxi_1) eg1 = eigen(K1) w = eg1$values
ww = w[-1]/w[-length(w)]
d = which(ww==min(ww[1:floor(0.75*p)]))

if (d > 1){
K1 = eg1$vectors[,1:d]%*%diag(eg1$values[1:d])%*%t(eg1$vectors[,1:d]); }else{ K1 = eg1$vectors[,1]%*%diag(eg1$values[1],1)%*%t(eg1$vectors[,1]);
}
K2 = S_yxi_1%*%t(S_yxi_2);
##(2) estimation of A and B
Geg = geigen::geigen(K2,K1);
evalues = Geg$values[which(Mod(Geg$values)<=10^5&Geg$values!=0)] Al = Geg$vectors[,which(Geg$values %in% evalues)] B = apply(t(S_yxi_1)%*%Al,2,l2s) Bl = t(MASS::ginv(B)) A = apply(S_yxi_1%*%Bl,2,l2s) } ##(3) estimation of Xt H = matrix(NA,p*q,d) for (ii in 1:d) { H[,ii] = B[,ii] %x% A[,ii] } f = Vec.tensor(Y)%*%H%*%MASS::ginv(t(H)%*%H) if(is.complex(A) == T || is.complex(B) == T ){ A = Complex2Real(A) B = Complex2Real(B) f = Complex2Real(f) } METHOD <- c("Estimation of matrix CP-factor model",paste("Method:",method)) names(d) <- "The estimated rank of the matrix CP-factor model" con = structure(list(A = A,B = B,f = f,Rank = d, method = METHOD), class = "mtscp") return(con) } if(method == "CP.Refined"){ ##(1) estimation of P,Q and d M1 = M2 = 0 dmax = round(min(p,q)*0.75) for (kk in 1:lag.k){ S_yxi_k = Autocov_xi_Y(Y,xi,lag.k = kk) M1 = M1 + S_yxi_k%*%t(S_yxi_k) M2 = M2 + t(S_yxi_k)%*%S_yxi_k } ev_M1 = eigen(M1) ev_M2 = eigen(M2) d1 = which.max(ev_M1$values[1:dmax]/ev_M1$values[2:(dmax+1)]) d2 = which.max(ev_M2$values[1:dmax]/ev_M2$values[2:(dmax+1)]) d = ifelse(p>q,d1,d2) if(!is.null(Rank)){ d = Rank } P = ev_M1$vectors[,1:d]
Q = ev_M2$vectors[,1:d] if(d == 1){ A = as.matrix(P) B = as.matrix(Q) f = vector() for (tt in 1:n) { f[tt] = t(A)%*%Y[tt,,]%*%B } f = as.matrix(f) }else{ ##(2) estimation of U and V Z = array(NA,dim = c(n,d,d)) for (tt in 1:n) { Z[tt,,] = t(P)%*%Y[tt,,]%*%Q } xi = est.xi(Z) S_Zxi_1 = Autocov_xi_Y(Z,xi,lag.k = 1) S_Zxi_2 = Autocov_xi_Y(Z,xi,lag.k = 2) vl = eigen(MASS::ginv(t(S_Zxi_1)%*%S_Zxi_1)%*%t(S_Zxi_1)%*%S_Zxi_2)$vectors ##MASS
ul = eigen(MASS::ginv(S_Zxi_1%*%t(S_Zxi_1))%*%S_Zxi_1%*%t(S_Zxi_2))$vectors U = apply(S_Zxi_1%*%vl,2,l2s) V = apply(t(S_Zxi_1)%*%ul,2,l2s) ##(3) estimation of A and B A = P%*%U B = Q%*%V ##(4) estimation of Xt W = matrix(NA,d^2,d) for (ii in 1:d) { W[,ii] = V[,ii]%x%U[,ii] } f = Vec.tensor(Z)%*%W%*%solve(t(W)%*%W) if(is.complex(A) == T || is.complex(B) == T ){ A = Complex2Real(A) B = Complex2Real(B) f = Complex2Real(f) } } METHOD <- c("Estimation of matrix CP-factor model",paste("Method:",method)) names(d) <- "The estimated rank of the matrix CP-factor model" con = structure(list(A = A,B = B,f = f,Rank = d, method = METHOD), class = "mtscp") return(con) } if(method == "CP.Unified"){ ##(1) estimation of P,Q and d1,d2 if(is.null(Rank)){ PQ_hat_tol = est.d1d2.PQ(Y,xi,K = lag.k) d1 = PQ_hat_tol$d1_hat
d2  = PQ_hat_tol$d2_hat d = NULL P = PQ_hat_tol$P_hat
Q   = PQ_hat_tol$Q_hat if(d1 == 1 || d2 == 1){d = d1*d2} }else{ d = Rank$d
d1  = Rank$d1 d2 = Rank$d2

PQ_hat_tol = est.PQ(Y,xi,d1,d2,K = lag.k)

P   = PQ_hat_tol$P_hat Q = PQ_hat_tol$Q_hat

}
##(2) estimation of W* = (v1*u1,v2*u2,...,vd*ud)H = WH
if(d1 == 1 & d2 == 1){
d = 1
f = vector()
for (tt in 1:n) {
f[tt] = t(P)%*%Y[tt,,]%*%Q
}
A = P
B = Q

# con = list(A = A,B = B,f = f,Rank = list(d=d,d1=d1,d2=d2))
METHOD <- c("Estimation of matrix CP-factor model",paste("Method:",method))
rank <- list(d=d,d1=d1,d2=d2)
names(rank) <- c("The estimated rank d of the matrix CP-factor model",
"The estimated rank d1 of the matrix CP-factor model",
"The estimated rankd d2 of the matrix CP-factor model")
con = structure(list(A = A,B = B,f = f,Rank = rank, method = METHOD),
class = "mtscp")

return(con)
}else{

if(is.null(d)){
W_hat_tol  =  est.d.Wf(Y,P,Q,Ktilde =  lag.ktilde)
d          =  W_hat_tol$d_hat W = W_hat_tol$W_hat
f          =  W_hat_tol$f_hat }else{ d = d W_hat_tol = est.Wf(Y,P,Q,d,Ktilde = lag.ktilde) W = W_hat_tol$W_hat
f          =  W_hat_tol$f_hat } ##(3) estimation of U and V if(d1 == 1 || d2 == 1){ Theta = NULL if(d1 == 1){ U = 1; V = W; warning("d1 equal to 1, V cannot be identified uniquely!") } if(d2 == 1){ U = W; V = 1; warning("d2 equal to 1, U cannot be identified uniquely!") } if(d1 == 1 & d2 == 1){ U = 1; V = 1; } U = as.matrix(U) V = as.matrix(V) }else{ UV_hat_tol = est.UV.JAD(W,d1,d2,d) U = UV_hat_tol$U
V          = UV_hat_tol$V Theta = UV_hat_tol$Theta
}

##(4) estimation of A and B
A = P%*%U
B = Q%*%V
METHOD <- c("Estimation of matrix CP-factor model",paste("Method:",method))
rank <- list(d=d,d1=d1,d2=d2)
names(rank) <- c("The estimated rank d of the matrix CP-factor model",
"The estimated rank d1 of the matrix CP-factor model",
"The estimated rankd d2 of the matrix CP-factor model")
con = structure(list(A = A,B = B,f = f,Rank = rank, method = METHOD),
class = "mtscp")

return(con)
}

}
}

rho2.loss = function(A_hat,A){
max(apply(1-(t(A_hat)%*%A)^2,2,min))
}

l2s = function(x){x/sqrt(sum(x^2))}

fnorm = function(x){sqrt(sum(x^2))}

Complex2Real = function(A){
REA = round(Re(A),8)
IMA = round(Im(A),8)

real.index  =  which(IMA[1,] == 0)

if(length(real.index) == 0){
complex_real  =  REA
complex_image =  IMA

complex_take  = which(duplicated(complex_real[1,]) == T)

real = as.matrix(complex_real[,complex_take])

img  = as.matrix(complex_image[,complex_take])

new_A = cbind(real,img)
}else{
real.vector   = REA[,real.index]

complex_real  =  REA[,-real.index]
complex_image =  IMA[,-real.index]

complex_take  = which(duplicated(complex_real[1,]) == T)

real = as.matrix(complex_real[,complex_take])

img  = as.matrix(complex_image[,complex_take])

new_A = cbind(real.vector,real,img)

colnames(new_A) = NULL
}

return(new_A)
}

Vec.tensor = function(Y){
p = dim(Y)[2];q = dim(Y)[3];

if(p == q & q == 1){
Y_tilde = apply(Y,1,FUN = as.vector)
}else{
Y_tilde = t(apply(Y,1,FUN = as.vector))
}
return(Y_tilde)

}

#' @title Data generate process of matrix CP-factor model
#' @description \code{DGP.CP()} function generate the matrix time series described in Chang et al. (2023):\deqn{{\bf{Y}}_t = {\bf A \bf X}_t{\bf B}^{'} +
#' {\boldsymbol{\epsilon}}_t, } where \eqn{{\bf X}_t = diag(x_{t,1},\ldots,x_{t,d})} is an \eqn{d \times d}
#' latent process, \eqn{{\bf A}} and \eqn{{\bf B}} are , respectively, \eqn{p
#' \times d} and \eqn{q \times d} unknown constant matrix, and \eqn{ {\boldsymbol{\epsilon}}_t }
#'  is a \eqn{p \times q} matrix white noise process.
#'
#' @param n  Integer. Sample size of \eqn{\bf Y_t}, \eqn{t=1,\ldots,n}.
#' @param p  Integer. Number of rows of \eqn{\bf Y_t}.
#' @param q  Integer. Number of columns of \eqn{\bf Y_t}.
#' @param d1  Integer. Rank of \eqn{\bf A}.
#' @param d2  Integer. Rank of \eqn{\bf B}.
#' @param d  Integer. Number of columns of \eqn{\bf A} and \eqn{\bf B}.
#'
#' @return A list containing the following
#'   components:
#'   \item{Y}{A \eqn{n \times p \times q} data array of \eqn{\bf Y_t}.}
#'   \item{S}{A \eqn{n \times p \times q} data array of \eqn{\bf S_t = \bf A \bf X_t \bf B'}.}
#'   \item{A}{A \eqn{p \times d} coefficient matrix.}
#'   \item{B}{A \eqn{q \times d} coefficient matrix.}
#'   \item{X}{A \eqn{n \times d \times d} data array of \eqn{\bf X_t}.}
#'   \item{P}{A \eqn{p \times d_1} orthogonal matrix such that \eqn{\bf A = \bf P \bf U}.}
#'   \item{Q}{A \eqn{q \times d_2} orthogonal matrix such that \eqn{\bf B = \bf Q \bf V}.}
#'   \item{U}{A \eqn{d_1 \times d} matrix such that \eqn{\bf A = \bf P \bf U}.}
#'   \item{V}{A \eqn{d_2 \times d} matrix such that \eqn{\bf B = \bf Q \bf V}.}
#'   \item{W}{A \eqn{d_1 d_2 \times d} matrix such that \eqn{\bf W = (\bf v_1 \otimes \bf u_1,\ldots,\bf v_d \otimes \bf u_d)}.}
#'   \item{Ws}{A \eqn{d_1 d_2 \times d} matrix. An orthogonal basis of \eqn{\bf W}.}
#'   \item{Xmat}{A \eqn{n \times d} data matrix of \eqn{diag(\bf X_t)}.}
#'   \item{Smat}{A \eqn{n \times pq} data matrix of \eqn{vec(\bf S_t)}.}
#' @references
#'   Chang, J., He, J., Yang, L. and Yao, Q.(2023). \emph{Modelling matrix time series via a tensor CP-decomposition}.
#'   Journal of the Royal Statistical Society Series B: Statistical Methodology, Vol. 85(1), pp.127--148.
#'
#' @examples
#' p = 10
#' q = 10
#' n = 400
#' d = d1 = d2 = 3
#' data <- DGP.CP(n,p,q,d1,d2,d)
#' Y = data$Y #' @export DGP.CP = function(n,p,q,d1,d2,d){ par_A = c(-3,3) par_B = c(-3,3) par_X = c(0.6,0.95) par_E = 1 Input = list(n = n,p = p,q = q,d1 = d1,d2 = d2,d = d) A_inl = matrix(runif(p*d,par_A[1],par_A[2]),p,d) B_inl = matrix(runif(q*d,par_B[1],par_B[2]),q,d) svd_A = svd(A_inl) P = svd_A$u[,1:d1]

A = (svd_A$u[,1:d1]) %*% diag(svd_A$d[1:d1],nrow = d1,ncol = d1) %*%t(svd_A$v[,1:d1]) U = t(P)%*%apply(A,2,l2s) A_s = P%*%U svd_B = svd(B_inl) Q = svd_B$u[,1:d2]
B = (svd_B$u[,1:d2]) %*% diag(svd_B$d[1:d2],nrow = d2,ncol = d2) %*% (t(svd_B$v[,1:d2])) V = t(Q)%*%apply(B,2,l2s) B_s = Q%*%V W = matrix(NA,d1*d2,d) for (ii in 1:d) { W[,ii] = V[,ii]%x%U[,ii] } X = array(0,c(n,d,d)) X_m = matrix(NA,n,d) signal = runif(d,-1,1) signal[signal>=0] <- 1 signal[signal< 0] <- -1 par_ar = runif(d,par_X[1],par_X[2]) for (ii in 1:d) { xd = arima.sim(model = list(ar = par_ar[ii]*signal[ii]),n = n) X_m[,ii] = xd*(apply(A, 2, fnorm)*apply(B, 2, fnorm))[ii] X[,ii,ii] <- xd*(apply(A, 2, fnorm)*apply(B, 2, fnorm))[ii] } S_m = X_m%*%t(W) W_star = svd(S_m)$v[,1:d]

Y = S = array(NA,dim = c(n,p,q))
for (tt in 1:n) {
S[tt,,] <- A_s%*%X[tt,,]%*%t(B_s)
Y[tt,,] <- S[tt,,] + matrix(rnorm(p*q,0,par_E),p,q)
}

return(list(Y      =  Y,
S      =  S,
A      =  A_s,
B      =  B_s,
X      =  X,
P      =  P,
Q      =  Q,
U      =  U,
V      =  V,
W      =  W,
Ws     =  W_star,
Xmat    =  X_m,
Smat    =  S_m,
Input  =  Input
))

}

Autocov_xi_Y = function(Y,xi,lag.k = k){

n = dim(Y)[1]
k = lag.k

Y_mean = 0
xi_mean = 0
for (ii in 1:n) {
Y_mean = Y_mean + Y[ii,,]
xi_mean = xi_mean + xi[ii]
}
Y_mean   = Y_mean/n
xi_mean = xi_mean/n

Sigma_Y_xi_k = 0
for (ii in (k+1):n) {
Sigma_Y_xi_k = Sigma_Y_xi_k + (Y[ii,,] - Y_mean)*(xi[ii-k] - xi_mean)
}
return(Sigma_Y_xi_k/(n-k))
}

est.xi  = function(Y, thresh_per = 0.99, d_max = 20){

n = dim(Y)[1];p = dim(Y)[2];q = dim(Y)[3];

xi.mat = Vec.tensor(Y)

if(n > p*q){
eig_xi.mat = eigen(MatMult(t(xi.mat),xi.mat))
cfr =  cumsum(eig_xi.mat$values)/sum(eig_xi.mat$values)
d_hat = min(which(cfr > thresh_per))
d_fin = min(d_max,d_hat)
w_hat = eig_xi.mat$vectors[,1:d_fin] xi.f = xi.mat%*%w_hat xi = rowMeans(xi.f) }else{ eig_xi.mat = eigen(MatMult(xi.mat,t(xi.mat))) cfr = cumsum(eig_xi.mat$values)/sum(eig_xi.mat$values) d_hat = min(which(cfr > thresh_per)) d_fin = min(d_max,d_hat) xi.f1 = as.matrix(eig_xi.mat$vectors[,1:d_fin])
weight = sqrt(eig_xi.mat$values[1:d_hat]) xi.f1 = xi.f1%*%diag(weight) xi = rowMeans(xi.f1) } return(xi) } est.d1d2.PQ = function(Y,xi,K = 10){ n = dim(Y)[1];p = dim(Y)[2];q = dim(Y)[3]; d2_list = d1_list =vector() M1 = M2 = 0 dmax = round(min(p,q)*0.75) P_list = Q_list = list() for (kk in 1:K){ S_yxi_k = Autocov_xi_Y(Y,xi,lag.k = kk) M1 = M1 + S_yxi_k%*%t(S_yxi_k) M2 = M2 + t(S_yxi_k)%*%S_yxi_k ev_M1 = eigen(M1) ev_M2 = eigen(M2) d1_list[kk] = which.max(ev_M1$values[1:dmax]/ev_M1$values[2:(dmax+1)]) d2_list[kk] = which.max(ev_M2$values[1:dmax]/ev_M2$values[2:(dmax+1)]) P_list[[kk]] = ev_M1$vectors[,1:(d1_list[kk])]
Q_list[[kk]] = ev_M2$vectors[,1:(d2_list[kk])] } d1_list[1] = 0 d2_list[1] = 0 d1_hat = d1_list[K] d2_hat = d2_list[K] P_hat = P_list[[K]] Q_hat = Q_list[[K]] return(list(d1_hat = d1_hat, d2_hat = d2_hat, P_hat = P_hat, Q_hat = Q_hat, d1_list = d1_list, d2_list = d2_list, P_list = P_list, Q_list = Q_list)) } est.PQ = function(Y,xi,d1,d2,K = 20){ n = dim(Y)[1];p = dim(Y)[2];q = dim(Y)[3]; M1 = M2 = 0 dmax = round(min(p,q)*0.75) P_list = Q_list = list() for (kk in 1:K){ S_yxi_k = Autocov_xi_Y(Y,xi,lag.k = kk) M1 = M1 + S_yxi_k%*%t(S_yxi_k) M2 = M2 + t(S_yxi_k)%*%S_yxi_k ev_M1 = eigen(M1) ev_M2 = eigen(M2) P_list[[kk]] = ev_M1$vectors[,1:(d1)]
Q_list[[kk]] = ev_M2$vectors[,1:(d2)] } P_hat = P_list[[K]] Q_hat = Q_list[[K]] return(list(P_hat = P_hat, Q_hat = Q_hat, P_list = P_list, Q_list = Q_list)) } est.d.Wf = function(Y,P,Q, Ktilde = 10){ n = dim(Y)[1];p = dim(Y)[2];q = dim(Y)[3]; d1 = NCOL(P);d2 = NCOL(Q); Z = array(NA,dim = c(n,d1,d2)) for (tt in 1:n) { Z[tt,,] = t(P)%*%Y[tt,,]%*%Q } Z_tilde = Vec.tensor(Z) M = 0 W_list = f_list = list() d_list = vector() dmax = d1*d2 dstar = max(d1,d2) for (kk in 1:Ktilde) { S_ztilde_k = sigmak(t(Z_tilde),as.matrix(colMeans(Z_tilde)),n = n, k= kk) M = M + S_ztilde_k%*%t(S_ztilde_k) ev_M = eigen(M) evalues = ev_M$values

d_list[kk] =  max(which.max(evalues[1:(dmax-1)]/evalues[2:(dmax)]),dstar)

W_list[[kk]] = ev_M$vectors[,1:(d_list[kk])] f_list[[kk]] = Z_tilde%*%(W_list[[kk]]) } d_list[1] = 0 d_hat = d_list[Ktilde] W_hat = W_list[[Ktilde]] f_hat = f_list[[Ktilde]] return(list(d_hat = d_hat, W_hat = W_hat, f_hat = f_hat, d_list = d_list, W_list = W_list, f_list = f_list)) } est.d.Wf.nPQ = function(Z, Ktilde = 10){ n = dim(Z)[1];d1 = dim(Z)[2];d2 = dim(Z)[3]; Z_tilde = Vec.tensor(Z) M = 0 W_list = f_list = list() d_list = vector() dmax = d1*d2 dstar = max(d1,d2) for (kk in 1:Ktilde){ S_ztilde_k = sigmak(t(Z_tilde),as.matrix(colMeans(Z_tilde)),n = n, k= kk) M = M + S_ztilde_k%*%t(S_ztilde_k) ev_M = eigen(M) evalues = ev_M$values

d_list[kk] =  max(which.max(evalues[1:(dmax-1)]/evalues[2:(dmax)]),dstar)

W_list[[kk]] = ev_M$vectors[,1:d_hat] f_list[[kk]] = Z_tilde%*%(W_list[[kk]]) } d_list[1] = 0 d_hat = d_list[Ktilde] W_hat = W_list[[Ktilde]] f_hat = f_list[[Ktilde]] return(list(d_hat = d_hat, W_hat = W_hat, f_hat = f_hat, d_list = d_list, W_list = W_list, f_list = f_list)) } est.Wf = function(Y,P,Q,d,Ktilde = 10){ n = dim(Y)[1];p = dim(Y)[2];q = dim(Y)[3]; d1 = NCOL(P);d2 = NCOL(Q); Z = array(NA,dim = c(n,d1,d2)) for (tt in 1:n) { Z[tt,,] = t(P)%*%Y[tt,,]%*%Q } Z_tilde = matrix(NA,n,d1*d2) for (tt in 1:n) { Z_tilde[tt,] = as.vector(Z[tt,,]) } M = 0 W_list = f_list = list() for (kk in 1:Ktilde){ S_ztilde_k = sigmak(t(Z_tilde),as.matrix(colMeans(Z_tilde)),n = n, k= kk) M = M + S_ztilde_k%*%t(S_ztilde_k) ev_M = eigen(M) W_list[[kk]] = ev_M$vectors[,1:d]
f_list[[kk]] = Z_tilde%*%(W_list[[kk]])
}

W_hat  =  W_list[[Ktilde]]
f_hat  =  f_list[[Ktilde]]
return(list(W_hat   = W_hat,
f_hat   = f_hat,
W_list  = W_list,
f_list  = f_list))
}

W_tilde_tol = array(NA,dim= c(d,d1,d2))

for (jj in 1:d){
W_tilde_i = matrix(NA,d1,d2)
for (pp in 1:d2){
for (mm in 1:d1) {
W_tilde_i[mm,pp] <- W[mm + (pp - 1)*d1,jj]
}
}
W_tilde_tol[jj,,] = W_tilde_i
}

P_tol = vector()
for (ss in 1:d) {
for(rr in ss:d){
if(ss == rr){
P_rs = minor_P(W_tilde_tol[rr,,],W_tilde_tol[ss,,],d1,d2)
}else{
P_rs = minor_P(W_tilde_tol[rr,,],W_tilde_tol[ss,,],d1,d2)
}
P_tol = cbind(P_tol,P_rs)
}
}

M_tol    = svd(P_tol)$v dt = NCOL(M_tol) M = M_tol[,(dt-d+1):dt] M_tensor = array(NA,dim = c(d,d,d)) eigen_gap = vector() for (ii in 1:d) { M_tensor[,,ii] = Vech2Mat_new(M[,ii],d) eigen_gap[ii] = min(abs(eigen(M_tensor[,,ii])$values))
}

##construct H_star
Ms = M_tensor[,,which.max(eigen_gap)]

HH0 = vector()
HH1 = vector()
HH2 = vector()

for (ii in 1:d) {
HH0 = cbind(HH0,c(M_tensor[,,ii]))
HH1 = cbind(HH1,c(solve(Ms)%*%M_tensor[,,ii]))
HH2 = cbind(HH2,c(M_tensor[,,ii]%*%solve(Ms)))
}

PP = (t(HH0)%*%HH2)%*%solve(t(HH1)%*%HH2 + t(HH2)%*%HH1)%*%(t(HH2)%*%HH0)

EVD = eigen(PP)
if( min(EVD$values) < 0){ R_NOJD = EVD$vectors
}else{
R_NOJD = sqrt(2)/2*(EVD$vectors)%*%diag((EVD$values)^{-1/2})%*%t(EVD$vectors) } M_1 = M%*%R_NOJD M_tensor_1 = array(NA,dim = c(d,d,d)) for (ii in 1:d) { M_tensor_1[,,ii] = Vech2Mat_new(M_1[,ii],d) } H = jointDiag::ffdiag(M_tensor_1)$B ## jointDiag::ffdiag

Theta = apply(MASS::ginv(H), 2, l2s)

Wt = W%*%Theta

U = matrix(NA,d1,d)
V = matrix(NA,d2,d)

for (jj in 1:d){
Wt_tilde_i = matrix(NA,d1,d2)
for (pp in 1:d2){
for (mm in 1:d1) {
Wt_tilde_i[mm,pp] <- Wt[mm + (pp - 1)*d1,jj]
}
}
svdi = svd(Wt_tilde_i)
U[,jj] = svdi$u[,1] V[,jj] = svdi$v[,1]
}

return(list(U = U,V = V,Theta = Theta))

}

est.UV.EVD = function(W,d1,d2,d){

W_tilde_tol = array(NA,dim= c(d,d1,d2))

for (jj in 1:d){
W_tilde_i = matrix(NA,d1,d2)
for (pp in 1:d2){
for (mm in 1:d1) {
W_tilde_i[mm,pp] <- W[mm + (pp - 1)*d1,jj]
}
}
W_tilde_tol[jj,,] = W_tilde_i
}

P_tol = vector()
for (ss in 1:d) {
for(rr in ss:d){
if(ss == rr){
P_rs = minor_P(W_tilde_tol[rr,,],W_tilde_tol[ss,,],d1,d2)
}else{
P_rs = minor_P(W_tilde_tol[rr,,],W_tilde_tol[ss,,],d1,d2)
}
P_tol = cbind(P_tol,P_rs)
}
}
Theta = 1 + 1i
M_tol    = svd(P_tol)$v dt = NCOL(M_tol) M_org = M_tol[,(dt-d+1):dt] M = M_org M_tensor = array(NA,dim = c(d,d,d)) for (ii in 1:d) { M_tensor[,,ii] = Vech2Mat_new(M[,ii],d) } L0 = L1 = 0 for (tt in 1:(d-1)) { L0 = L0 + M_tensor[,,tt] L1 = L1 + M_tensor[,,tt + 1] } L0 = L0/(d-1) L1 = L1/(d-1) theta.l = eigen(solve(L1)%*%L0)$vectors

Theta = apply(L0%*%theta.l,2,l2s)

if(is.complex(Theta)){
Theta = Complex2Real(Theta)
}

Wt = W%*%Theta

U = matrix(NA,d1,d)
V = matrix(NA,d2,d)

for (jj in 1:d){
Wt_tilde_i = matrix(NA,d1,d2)
for (pp in 1:d2){
for (mm in 1:d1) {
Wt_tilde_i[mm,pp] <- Wt[mm + (pp - 1)*d1,jj]
}
}
svdi = svd(Wt_tilde_i)
U[,jj] = svdi$u[,1] V[,jj] = svdi$v[,1]
}

return(list(U = U,V = V,Theta = Theta))

}


Try the HDTSA package in your browser

Any scripts or data that you put into this service are public.

HDTSA documentation built on Sept. 11, 2024, 5:49 p.m.