# R/tenAR.R In tensorTS: Factor and Autoregressive Models for Tensor Time Series

#### Documented in matAR.RR.estmatAR.RR.semplotmplot.acftenAR.esttenAR.predicttenAR.sim

###Functions of Autoregressive Models

#' Generate TenAR(p) tensor time series
#'
#' Simulate from the TenAR(p) model.
#'@name tenAR.sim
#'@rdname tenAR.sim
#'@aliases tenAR.sim
#'@export
#'@import tensor rTensor
#'@importFrom MASS ginv
#'@importFrom stats rnorm
#'@importFrom pracma randortho
#'@importFrom Matrix rankMatrix
#'@param t length of output series, a strictly positive integer.
#'@param dim dimension of the tensor at each time.
#'@param R Kronecker rank for each lag.
#'@param P autoregressive order.
#'@param rho spectral radius of coefficient matrix \eqn{\Phi}.
#'@param cov covariance matrix of the error term: diagonal ("iid"), separable ("mle"), random ("svd").
#'@param A coefficient matrices. If not provided, they are randomly generated according to given \code{dim}, \code{R}, \code{P} and \code{rho}.
#' It is a multi-layer list, the first layer for the lag \eqn{1 \le i \le P}, the second the term \eqn{1 \le r \le R}, and the third the mode \eqn{1 \le k \le K}.
#'@param Sig only if \code{cov=mle}, a list of initial values of \eqn{\Sigma_1,\ldots,\Sigma_K}. The default are identity matrices.
#'@return A tensor-valued time series generated by the TenAR(p) model.
#'@examples
#' set.seed(123)
#' dim <- c(3,3,3)
#' xx <- tenAR.sim(t=500, dim, R=2, P=1, rho=0.5, cov='iid')
tenAR.sim <- function(t, dim, R, P, rho, cov, A=NULL, Sig=NULL){
if (is.null(A)){A <- tenAR.A(dim, R, P, rho)}
K <- length(A[[1]][[1]])
dim <- c()
for (i in c(1:K)){
dim[i] <- nrow(A[[1]][[1]][[i]])
}
x <- array(0, c(t+500,prod(dim)))
for (l in c(1:P)){
x[l,] <- rnorm(prod(dim))
}

if (cov == "mle"){
if (is.null(Sig)){
Sig.true <- lapply(1:K, function(i){
Q <- pracma::randortho(dim[i])
D <- abs(diag(rnorm(dim[i])))
Q %*% D %*% t(Q)
})
Sig.true <- fro.rescale(list(Sig.true))[[1]]
Sig.true.sqrt <- lapply(1:K, function(i){
sqrtm(Sig.true[[i]])$B }) } else { Sig.true = Sig Sig.true.sqrt <- lapply(1:K, function(i){ sqrtm(Sig.true[[i]])$B
})
}
}

if (cov == "svd"){
Q <- pracma::randortho(prod(dim))
D <- sqrt(abs(diag(rnorm(prod(dim)))))
E <- Q %*% D
}

for (i in c((P+1):(500 + t))){ # burning number = 500
if (cov == "iid"){
e <- rnorm(prod(dim), mean=0, sd=1)
} else if (cov == "mle"){
e <- kronecker_list(rev(Sig.true.sqrt)) %*% rnorm(prod(dim))
} else if (cov == "svd"){
e <- E %*% rnorm(prod(dim))
} else {
}
temp = 0
for (l in c(1:P)){
if (R[l] == 0) next
phi <-  Reduce("+", lapply(1:R[l], function(j) {rTensor::kronecker_list(rev(A[[l]][[j]]))}))
temp = temp + phi %*% x[i-l, ]
}
x[i,] <-  temp + e
}
return(array(x[501:(500+t),], c(t, dim)))
}

#' Estimation for Autoregressive Model of Tensor-Valued Time Series
#'
#' Estimation function for tensor autoregressive models. Methods include
#' projection (PROJ), Least Squares (LSE), maximum likelihood estimation (MLE)
#' and vector autoregressive model (VAR), as determined by the value of \code{method}.
#'@details
#' Tensor autoregressive model (of autoregressive order one) has the form:
#' \deqn{X_t = \sum_{r=1}^R X_{t-1} \times_{1}  A_1^{(r)} \times_{2}  \cdots \times_{K} A_K^{(r)} + E_t,}
#' where \eqn{A_k^{(r)}} are \eqn{d_k \times d_k} coefficient matrices, \eqn{k=1,\cdots,K}, and \eqn{E_t} is a tensor white noise. \eqn{R} is the Kronecker rank.
#' The model of autoregressive order \eqn{P} takes the form
#' \deqn{X_t = \sum_{i=1}^{P} \sum_{r=1}^{R_i} X_{t-i} \times_{1} A_{1}^{(ir)} \times_{2}  \cdots \times_{K} A_{K}^{(ir)} + E_t.}
#' For the "MLE" method, we also assume,
#'  \deqn{\mathrm{Cov}(\mathrm{vec}(E_t))= \Sigma_K \otimes \Sigma_{K-1} \otimes \cdots \otimes \Sigma_1,}
#'@name tenAR.est
#'@rdname tenAR.est
#'@aliases tenAR.est
#'@usage tenAR.est(xx,R=1,P=1,method="LSE",init.A=NULL,init.sig=NULL,niter=150,tol=1e-6)
#'@export
#'@import tensor rTensor
#'@importFrom methods new
#'@param xx \eqn{T \times d_1 \times \cdots \times d_K} tensor-valued time series, \eqn{T} is the length of the series.
#'@param method character string, specifying the type of the estimation method to be used. \describe{
#'  \item{\code{"PROJ",}}{Projection method.}
#'  \item{\code{"LSE",}}{Least squares.}
#'  \item{\code{"MLE",}}{MLE under a separable cov(vec(\eqn{E_t})).}
#'  \item{\code{"VAR",}}{VAR(\eqn{P}) model for the \eqn{\mathrm{vec}(E_t)}.}
#'}
#'@param R Kronecker rank for each lag.
#'@param P Autoregressive order.
#'@param init.A initial values of coefficient matrices \eqn{A_k^{(ir)}} in estimation algorithms, which is a multi-layer list such that
#' the first layer for the lag \eqn{1 \le i \le P}, the second the term \eqn{1 \le r \le R}, and the third the mode \eqn{1 \le k \le K}.
#' See "Details". By default, we use PROJ estimators as initial values.
#'@param init.sig only if \code{method=MLE}, a list of initial values of \eqn{\Sigma_1,\ldots,\Sigma_K}. The default are identity matrices.
#'@param niter maximum number of iterations if error stays above \code{tol}.
#'@param tol error tolerance in terms of the Frobenius norm.
#'@return return a list containing the following:\describe{
#'\item{\code{A}}{a list of estimated coefficient matrices \eqn{A_k^{(ir)}}. It is a multi-layer list,
#' the first layer for the lag \eqn{1 \le i \le P}, the second the term \eqn{1 \le r \le R}, and the third the mode \eqn{1 \le k \le K}. See "Details".}
#'\item{\code{SIGMA}}{only if \code{method=MLE}, a list of estimated \eqn{\Sigma_1,\ldots,\Sigma_K}.}
#'\item{\code{res}}{residuals}
#'\item{\code{Sig}}{covariance matrix \eqn{\mathrm{cov}(\mathrm{vec}({E_t}))}.}
#'\item{\code{cov}}{grand covariance matrix of all estimated entries of \eqn{A_k^{(ir)}}}
#'\item{\code{sd}}{standard errors of the coefficient matrices \eqn{A_k^{(ir)}}, returned as a list aligned with \code{A}.}
#'\item{\code{niter}}{number of iterations.}
#'\item{\code{BIC}}{value of extended Bayesian information criterion.}
#'}
#'@references
#'Rong Chen, Han Xiao, and Dan Yang. "Autoregressive models for matrix-valued time series". Journal of Econometrics, 2020.
#'@examples
#' set.seed(333)
#'
#' # case 1: tensor-valued time series
#'
#' dim <- c(2,2,2)
#' xx <- tenAR.sim(t=100, dim, R=2, P=1, rho=0.5, cov='iid')
#' est <- tenAR.est(xx, R=2, P=1, method="LSE")
#' A <- est$A # A is a multi-layer list #' #' length(A) == 1 # TRUE, since the order P = 1 #' length(A[[1]]) == 2 # TRUE, since the number of terms R = 2 #' length(A[[1]][[1]]) == 3 # TRUE, since the mode K = 3 #' #' #' # case 2: matrix-valued time series #' #' dim <- c(2,2) #' xx <- tenAR.sim(t=100, dim, R=2, P=1, rho=0.5, cov='iid') #' est <- tenAR.est(xx, R=2, P=1, method="LSE") #' A <- est$A # A is a multi-layer list
#'
#' length(A) == 1 # TRUE, since the order P = 1
#' length(A[[1]]) == 2 # TRUE, since the number of terms R = 2
#' length(A[[1]][[1]]) == 2 # TRUE, since the mode K = 2
tenAR.est <- function(xx, R=1, P=1, method="LSE", init.A=NULL, init.sig=NULL, niter=150, tol=1e-6){
if (identical("PROJ", method)) {
tenAR.PROJ(xx, R, P)
} else if (identical("LSE", method)) {
tenAR.LS(xx, R, P, init.A, niter, tol, print.true=FALSE)
} else if (identical("MLE", method)) {
tenAR.MLE(xx, R, P, init.A, init.sig, niter, tol, print.true=FALSE)
} else if (identical("VAR", method)) {
tenAR.VAR(xx, P)
} else {
stop("Please specify the type you want to use. See manuals or run ?tenAR for details.")
}
}

#' Estimation for Reduced Rank MAR(1) Model
#'
#' Estimation of the reduced rank MAR(1) model, using least squares (RRLSE) or MLE (RRMLE), as determined by the value of \code{method}.
#'@details
#' The reduced rank MAR(1) model takes the form:
#' \deqn{X_t =  A_1 X_{t-1} A_2^{\prime} + E_t,}
#' where \eqn{A_i} are \eqn{d_i \times d_i} coefficient matrices of ranks \eqn{\mathrm{rank}(A_i) = k_i \le d_i}, \eqn{i=1,2}. For the MLE method we also assume
#'  \deqn{\mathrm{Cov}(\mathrm{vec}(E_t))=\Sigma_2 \otimes \Sigma_1}
#'
#'@name matAR.RR.est
#'@rdname matAR.RR.est
#'@aliases matAR.RR.est
#'@usage matAR.RR.est(xx, method, A1.init=NULL, A2.init=NULL,Sig1.init=NULL,Sig2.init=NULL,
#'k1=NULL, k2=NULL, niter=100,tol=1e-6)
#'@export
#'@import tensor rTensor expm
#'@importFrom MASS ginv
#'@param xx \eqn{T \times d_1 \times d_2} matrix-valued time series, \eqn{T} is the length of the series.
#'@param method character string, specifying the method of the estimation to be used. \describe{
#'  \item{\code{"RRLSE",}}{Least squares.}
#'  \item{\code{"RRMLE",}}{MLE under a separable cov(vec(\eqn{E_t})).}
#'}
#'@param A1.init initial value of \eqn{A_1}. The default is the identity matrix.
#'@param A2.init  initial value of \eqn{A_2}. The default is the identity matrix.
#'@param Sig1.init only if \code{method=RRMLE}, initial value of \eqn{\Sigma_1}. The default is the identity matrix.
#'@param Sig2.init only if \code{method=RRMLE}, initial value of \eqn{\Sigma_2}. The default is the identity matrix.
#'@param k1  rank of \eqn{A_1}, a positive integer.
#'@param k2  rank of \eqn{A_2}, a positive integer.
#'@param niter maximum number of iterations if error stays above \code{tol}.
#'@param tol tolerance in terms of the Frobenius norm.
#'@param niter maximum number of iterations if error stays above \code{tol}.
#'@param tol relative Frobenius norm error tolerance.
#'@return return a list containing the following:\describe{
#'\item{\code{A1}}{estimator of \eqn{A_1}, a \eqn{d_1} by \eqn{d_1} matrix}
#'\item{\code{A2}}{estimator of \eqn{A_2}, a \eqn{d_2} by \eqn{d_2} matrix}
#'\item{\code{Sig1}}{only if \code{method=MLE}, when \eqn{\mathrm{Cov}(\mathrm{vec}(E_t))=\Sigma_2 \otimes \Sigma_1}.}
#'\item{\code{Sig2}}{only if \code{method=MLE}, when \eqn{\mathrm{Cov}(\mathrm{vec}(E_t))=\Sigma_2 \otimes \Sigma_1}.}
#'\item{\code{res}}{residuals.}
#'\item{\code{Sig}}{covariance matrix cov(vec(\eqn{E_t}))}
#'\item{\code{cov}}{covariance matrix \eqn{\hat A_1} and \eqn{\hat A_2}. If \code{method=RRLSE} or \code{method=RRMLE}, then it is a list containing \describe{
#'  \item{\code{Sigma}}{asymptotic covariance matrix of (vec( \eqn{\hat A_1}),vec(\eqn{\hat A_2^T})).}
#'  \item{\code{Theta1.u}, \code{Theta1.v}}{asymptotic covariance matrix of vec(\eqn{\hat U_1}), vec(\eqn{\hat V_1}).}
#'  \item{\code{Theta2.u}, \code{Theta2.v}}{asymptotic covariance matrix of vec(\eqn{\hat U_2}), vec(\eqn{\hat V_2}).}
#'}}
#'\item{\code{sd}}{standard errors of \eqn{\hat A_1} and \eqn{\hat A_2}, returned in a list aligned with \eqn{\hat A_1} and \eqn{\hat A_2}.}
#'\item{\code{niter}}{number of iterations.}
#'\item{\code{BIC}}{value of the extended Bayesian information criterion.}
#'}
#'@references
#'Reduced Rank Autoregressive Models for Matrix Time Series, by Han Xiao, Yuefeng Han, Rong Chen and Chengcheng Liu.
#'@examples
#' set.seed(333)
#' dim <- c(3,3)
#' xx <- tenAR.sim(t=500, dim, R=2, P=1, rho=0.5, cov='iid')
#' est <- matAR.RR.est(xx, method="RRLSE", k1=1, k2=1)
matAR.RR.est <- function(xx, method, A1.init=NULL, A2.init=NULL, Sig1.init=NULL, Sig2.init=NULL,k1=NULL, k2=NULL, niter=100,tol=1e-6){
if (identical("PROJ", method)) {
MAR1.PROJ(xx) # just keep it there
} else if (identical("LSE", method)) {
MAR1.LS(xx) # just keep it there
} else if (identical("MLE", method)) {
MAR1.MLE(xx) # just keep it there
} else if (identical("VAR", method)) {
tenAR.VAR(xx, P=1)
} else if (identical("RRLSE", method)) {
MAR1.RR(xx, k1, k2, niter, tol, A1.init, A2.init)
} else if (identical("RRMLE", method)) {
MAR1.CC(xx, k1, k2, A1.init, A2.init, Sigl.init=Sig1.init, Sigr.init=Sig2.init, niter, tol)
} else {
stop("Please specify the method in MAR.")
}
}

#' Asymptotic Covariance Matrix of One-Term Reduced rank MAR(1) Model
#'
#' Asymptotic covariance matrix of the reduced rank MAR(1) model. If \code{Sigma1} and \code{Sigma2} is provided in input,
#' we assume a separable covariance matrix, Cov(vec(\eqn{E_t})) = \eqn{\Sigma_2 \otimes \Sigma_1}.
#'@name matAR.RR.se
#'@rdname matAR.RR.se
#'@aliases matAR.RR.se
#'@usage matAR.RR.se(A1,A2,k1,k2,method,Sigma.e=NULL,Sigma1=NULL,Sigma2=NULL,RU1=diag(k1),
#'RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100)
#'@import tensor rTensor expm
#'@importFrom MASS ginv
#'@export
#'@param A1 left coefficient matrix.
#'@param A2 right coefficient matrix.
#'@param k1 rank of \eqn{A_1}.
#'@param k2 rank of \eqn{A_2}.
#'@param method character string, specifying the method of the estimation to be used. \describe{
#'  \item{\code{"RRLSE",}}{Least squares.}
#'  \item{\code{"RRMLE",}}{MLE under a separable cov(vec(\eqn{E_t})).}
#'}
#'@param Sigma.e only if \code{method} = "RRLSE". Cov(vec(\eqn{E_t})) = Sigma.e: covariance matrix of dimension \eqn{(d_1 d_2) \times (d_1 d_2)}
#'@param Sigma1,Sigma2 only if \code{method} = "RRMLE". Cov(vec(\eqn{E_t})) = \eqn{\Sigma_2 \otimes \Sigma_1}. \eqn{\Sigma_i} is \eqn{d_i \times d_i}, \eqn{i=1,2}.
#'@param RU1,RV1,RU2,RV2 orthogonal rotations of \eqn{U_1,V_1,U_2,V_2}, e.g., new_U1=U1 RU1
#'@param mpower truncate the VMA(\eqn{\infty}) representation of vec(\eqn{X_t}) at \code{mpower} for the purpose of calculating the autocovariances. The default is 100.
#'@return a list containing the following:\describe{
#'\item{\code{Sigma}}{asymptotic covariance matrix of (vec(\eqn{\hat A_1}),vec(\eqn{\hat A_2^T})).}
#'\item{\code{Theta1.u}}{asymptotic covariance matrix of vec(\eqn{\hat U_1}).}
#'\item{\code{Theta1.v}}{asymptotic covariance matrix of vec(\eqn{\hat V_1}).}
#'\item{\code{Theta2.u}}{asymptotic covariance matrix of vec(\eqn{\hat U_2}).}
#'\item{\code{Theta2.v}}{asymptotic covariance matrix of vec(\eqn{\hat V_2}).}
#'}
#'@references
#'Han Xiao, Yuefeng Han, Rong Chen and Chengcheng Liu, Reduced Rank Autoregressive Models for Matrix Time Series.
matAR.RR.se <- function(A1,A2,k1,k2,method,Sigma.e=NULL,Sigma1=NULL,Sigma2=NULL,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100){
if (method == "RRLSE"){
return(MAR1.RRLS.SE(A1,A2,k1,k2,Sigma.e,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100))
} else if (method == "RRMLE") {
return(MAR1.RRCC.SE(A1,A2,k1,k2,Sigma1,Sigma2,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100))
} else {
}
}

MAR1.RRLS.SE <- function(A1,A2,k1,k2,Sigma.e,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100){
# iterative least square
# X_t = A1 X_{t-1} A2^T + E_t
# RU1,RV1,RU2,RV2: rotation of U1,V1,U2,V2, e.g., new_U1=U1 RU1
# Cov(vec(E_t)) = Sigma.e : of dimension d1d2\times d1d2
# truncate vec(X_t) to mpower term, i.e. VMA(mpower)
# return Sigma.LS: asymptotic covariance matrix of (vec(\hat A1),vec(\hat A2^T))
#  Theta1.LS.u, Theta1.LS.v: asymptotic covariance matrix of vec(\hat U1), vec(\hat V1)
#  Theta2.LS.u, Theta2.LS.v: asymptotic covariance matrix of vec(\hat U2), vec(\hat V2)
d1 <- dim(A1)[1]
d2 <- dim(A2)[1]
Sigma.ee <- array(Sigma.e,c(d1,d2,d1,d2))
B <- kronecker(A2,A1)
#Sigma.x.vec, Sigma.xx: covariance matrix of vec(X_t)
Sigma.x.vec <- Sigma.e
for(i in 1:mpower){
Sigma.x.vec <- Sigma.x.vec+(B%^%i)%*%Sigma.e%*%t(B%^%i)
}
Sigma.xx <- array(Sigma.x.vec,c(d1,d2,d1,d2))
Gamma1 <- tensor(Sigma.xx,t(A1)%*%A1,c(1,3),c(1,2))
Gamma2 <- tensor(Sigma.xx,t(A2)%*%A2,c(2,4),c(1,2))
Gamma3 <- kronecker(diag(1,d1),A1)%*%matrix(aperm(Sigma.xx,c(3,1,4,2)),d1^2,d2^2)%*% kronecker(t(A2),diag(1,d2))
H <- matrix(NA,d1^2+d2^2,d1^2+d2^2)
H[1:d1^2,1:d1^2] <- as.vector(A1)%*%t(as.vector(A1)) + kronecker(Gamma2,diag(1,d1))
H[1:d1^2,(d1^2+1):(d1^2+d2^2)] <- Gamma3
H[(d1^2+1):(d1^2+d2^2),1:d1^2] <- t(Gamma3)
H[(d1^2+1):(d1^2+d2^2),(d1^2+1):(d1^2+d2^2)] <- kronecker(diag(1,d2),Gamma1)
H.inv <- solve(H)
D1 <- Gamma2%*%t(A1)%*% ginv(A1%*%Gamma2%*%t(A1)) %*%A1
D2 <- Gamma1%*%t(A2)%*% ginv(A2%*%Gamma1%*%t(A2)) %*%A2
P1 <- A1%*%ginv(t(A1)%*%A1)%*%t(A1)
P2 <- A2%*%ginv(t(A2)%*%A2)%*%t(A2)
Sigma.Q <- matrix(NA,d1^2+d2^2,d1^2+d2^2)
M1 <- kronecker(t(A2),P1)%*%Sigma.e%*%kronecker(A2,P1)
M2 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e%*%kronecker(A2,P1)
M3 <- kronecker(t(A2),P1)%*%Sigma.e%*%kronecker(A2,diag(d1)-P1)
M4 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e%*%kronecker(A2,diag(d1)-P1)
Sigma.Q1 <- Sigma.Q2 <- Sigma.Q3 <- Sigma.Q4 <- matrix(0,d1^2,d1^2)
for(i in 1:d1){
for(j in 1:d1){
for(k in 1:d1){
for(l in 1:d1){
Sigma.Q1[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M1[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
Sigma.Q2[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M2[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
Sigma.Q3[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M3[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
Sigma.Q4[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M4[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
}
}
}
}
Sigma.Q[1:d1^2,1:d1^2] <- Sigma.Q1+kronecker(D1,diag(d1))%*%Sigma.Q2+Sigma.Q3%*%kronecker(t(D1),diag(d1))+
kronecker(D1,diag(d1))%*%Sigma.Q4%*%kronecker(t(D1),diag(d1))
M1 <- kronecker(P2,t(A1))%*%Sigma.e%*%kronecker(P2,A1)
M2 <- kronecker(P2,t(A1))%*%Sigma.e%*%kronecker(diag(d2)-P2,A1)
M3 <- kronecker(diag(d2)-P2,t(A1))%*%Sigma.e%*%kronecker(P2,A1)
M4 <- kronecker(diag(d2)-P2,t(A1))%*%Sigma.e%*%kronecker(diag(d2)-P2,A1)
Sigma.Q1 <- Sigma.Q2 <- Sigma.Q3 <- Sigma.Q4 <- matrix(0,d2^2,d2^2)
for(i in 1:d2){
for(j in 1:d2){
for(k in 1:d2){
for(l in 1:d2){
Sigma.Q1[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M1[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
Sigma.Q2[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M2[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
Sigma.Q3[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M3[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
Sigma.Q4[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M4[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
}
}
}
}
Sigma.Q[(d1^2+1):(d1^2+d2^2),(d1^2+1):(d1^2+d2^2)] <- Sigma.Q1+Sigma.Q2%*%kronecker(diag(d2),t(D2))+
kronecker(diag(d2),D2)%*%Sigma.Q3+kronecker(diag(d2),D2)%*%Sigma.Q4%*%kronecker(diag(d2),t(D2))
M1 <- kronecker(t(A2),P1)%*%Sigma.e%*%kronecker(P2,A1)
M2 <- kronecker(t(A2),P1)%*%Sigma.e%*%kronecker(diag(d2)-P2,A1)
M3 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e%*%kronecker(P2,A1)
M4 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e%*%kronecker(diag(d2)-P2,A1)
Sigma.Q1 <- Sigma.Q2 <- Sigma.Q3 <- Sigma.Q4 <- matrix(0,d1^2,d2^2)
for(i in 1:d1){
for(j in 1:d1){
for(k in 1:d2){
for(l in 1:d2){
Sigma.Q1[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M1[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
Sigma.Q2[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M2[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
Sigma.Q3[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M3[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
Sigma.Q4[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M4[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
}
}
}
}
Sigma.Q[1:d1^2,(d1^2+1):(d1^2+d2^2)] <- Sigma.Q1+Sigma.Q2%*%kronecker(diag(d2),t(D2))+
kronecker(D1,diag(d1))%*%Sigma.Q3+kronecker(D1,diag(d1))%*%Sigma.Q4%*%kronecker(diag(d2),t(D2))
Sigma.Q[(d1^2+1):(d1^2+d2^2),1:d1^2] <- t(Sigma.Q[1:d1^2,(d1^2+1):(d1^2+d2^2)])
Sigma.LS <- H.inv%*%Sigma.Q%*%H.inv
Sigma.LS11 <- Sigma.LS[1:d1^2,1:d1^2]
Sigma.LS11.t <- array(Sigma.LS11,c(d1,d1,d1,d1))
Sigma.LS11.t <- aperm(Sigma.LS11.t,c(2,1,4,3))
Sigma.LS11.t <- matrix(Sigma.LS11.t, d1^2, d1^2)
Sigma.LS22.t <- Sigma.LS[(d1^2+1):(d1^2+d2^2),(d1^2+1):(d1^2+d2^2)]
Sigma.LS22 <- array(Sigma.LS22.t,c(d2,d2,d2,d2))
Sigma.LS22 <- aperm(Sigma.LS22,c(2,1,4,3))
Sigma.LS22 <- matrix(Sigma.LS22, d2^2, d2^2)
svd.A1 <- svd(A1)
#k1 <- length(svd.A1$d[svd.A1$d>1e-10])
D1 <- diag(c(svd.A1$d[1:k1],1))[1:k1,1:k1] U1 <- svd.A1$u[,1:k1]
V1 <- svd.A1$v[,1:k1] if(k1<d1){ U1c <- svd.A1$u[,(k1+1):d1]
V1c <- svd.A1$v[,(k1+1):d1] }else if(k1==d1){ U1c <- 0 V1c <- 0 } #U1c <- svd.A1$u[,(k1+1):d1]  # original version : didn't consider if k1=d1
#V1c <- svd.A1$v[,(k1+1):d1] e1 <- diag(d1) J1 <- matrix(0,d1^2,d1^2) for(i in 1:d1){ J1[,((i-1)*d1+1):(i*d1)] <- kronecker(diag(d1),e1[,i]) } C1 <- kronecker(A1,diag(d1)) + kronecker(diag(d1),A1)%*%J1 Sigma.LS.1u <- C1%*%Sigma.LS11%*%t(C1) C1.t <- kronecker(t(A1),diag(d1)) + kronecker(diag(d1),t(A1))%*%J1 Sigma.LS.1v <- C1.t%*%Sigma.LS11.t%*%t(C1.t) e1 <- diag(k1) L1 <- matrix(0,k1^2,k1) for(i in 1:k1){ L1[,i] <- kronecker(e1[,i],e1[,i]) } if(k1<d1){ R1u <- cbind(kronecker(diag(k1),U1), kronecker(diag(k1),U1c)) %*% rbind(solve(kronecker(D1^2,diag(k1)) - kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(U1),t(U1)), kronecker(diag(1/c(svd.A1$d[1:k1],1))[1:k1,1:k1]^2%*%t(U1), t(U1c))  )
R1v <- cbind(kronecker(diag(k1),V1), kronecker(diag(k1),V1c)) %*% rbind(solve(kronecker(D1^2,diag(k1)) -
kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(V1),t(V1)),
kronecker(diag(1/c(svd.A1$d[1:k1],1))[1:k1,1:k1]^2%*%t(V1), t(V1c)) ) }else if(k1==d1){ R1u <- kronecker(diag(k1),U1) %*% solve(kronecker(D1^2,diag(k1)) - kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(U1),t(U1)) R1v <- kronecker(diag(k1),V1) %*% solve(kronecker(D1^2,diag(k1)) - kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(V1),t(V1)) } # R1u <- cbind(kronecker(diag(k1),U1), kronecker(diag(k1),U1c)) %*% rbind(solve(kronecker(D1^2,diag(k1)) - # kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(U1),t(U1)), # kronecker(diag(1/c(svd.A1$d[1:k1],1))[1:k1,1:k1]^2%*%t(U1), t(U1c))  )
# R1v <- cbind(kronecker(diag(k1),V1), kronecker(diag(k1),V1c)) %*% rbind(solve(kronecker(D1^2,diag(k1)) -
#                                                                                 kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(V1),t(V1)),
#                                                                         kronecker(diag(1/c(svd.A1$d[1:k1],1))[1:k1,1:k1]^2%*%t(V1), t(V1c)) ) Theta1.LS.u <- R1u%*% Sigma.LS.1u %*%t(R1u) Theta1.LS.v <- R1v%*% Sigma.LS.1v %*%t(R1v) Theta1.LS.u <- kronecker(t(RU1),diag(d1))%*%Theta1.LS.u Theta1.LS.v <- kronecker(t(RV1),diag(d1))%*%Theta1.LS.v svd.A2 <- svd(A2) #k2 <- length(svd.A2$d[svd.A2$d>1e-10]) D2 <- diag(c(svd.A2$d[1:k2],1))[1:k2,1:k2]
U2 <- svd.A2$u[,1:k2] V2 <- svd.A2$v[,1:k2]
if(k2<d2){
U2c <- svd.A2$u[,(k2+1):d2] V2c <- svd.A2$v[,(k2+1):d2]
}else if(k2==d2){
U2c <- 0
V2c <- 0
}
#U2c <- svd.A2$u[,(k2+1):d2] #V2c <- svd.A2$v[,(k2+1):d2]
e2 <- diag(d2)
J2 <- matrix(0,d2^2,d2^2)
for(i in 1:d2){
J2[,((i-1)*d2+1):(i*d2)] <- kronecker(diag(d2),e2[,i])
}
C2 <- kronecker(A2,diag(d2)) + kronecker(diag(d2),A2)%*%J2
Sigma.LS.2u <- C2%*%Sigma.LS22%*%t(C2)
C2.t <- kronecker(t(A2),diag(d2)) + kronecker(diag(d2),t(A2))%*%J2
Sigma.LS.2v <- C2.t%*%Sigma.LS22.t%*%t(C2.t)
e2 <- diag(k2)
L2 <- matrix(0,k2^2,k2)
for(i in 1:k2){
L2[,i] <- kronecker(e2[,i],e2[,i])
}

if(k2<d2){
R2u <- cbind(kronecker(diag(k2),U2), kronecker(diag(k2),U2c)) %*% rbind(solve(kronecker(D2^2,diag(k2)) -
kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(U2),t(U2)),
kronecker(diag(1/c(svd.A2$d[1:k2],1))[1:k2,1:k2]^2%*%t(U2), t(U2c)) ) R2v <- cbind(kronecker(diag(k2),V2), kronecker(diag(k2),V2c)) %*% rbind(solve(kronecker(D2^2,diag(k2)) - kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(V2),t(V2)), kronecker(diag(1/c(svd.A2$d[1:k2],1))[1:k2,1:k2]^2%*%t(V2), t(V2c))  )
}else if(k2==d2){
R2u <- kronecker(diag(k2),U2) %*% solve(kronecker(D2^2,diag(k2)) -
kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(U2),t(U2))
R2v <- kronecker(diag(k2),V2) %*% solve(kronecker(D2^2,diag(k2)) -
kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(V2),t(V2))

}

# R2u <- cbind(kronecker(diag(k2),U2), kronecker(diag(k2),U2c)) %*% rbind(solve(kronecker(D2^2,diag(k2)) -
#                                                                                 kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(U2),t(U2)),
#                                                                         kronecker(diag(1/c(svd.A2$d[1:k2],1))[1:k2,1:k2]^2%*%t(U2), t(U2c)) ) # R2v <- cbind(kronecker(diag(k2),V2), kronecker(diag(k2),V2c)) %*% rbind(solve(kronecker(D2^2,diag(k2)) - # kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(V2),t(V2)), # kronecker(diag(1/c(svd.A2$d[1:k2],1))[1:k2,1:k2]^2%*%t(V2), t(V2c))  )
Theta2.LS.u <- R2u%*% Sigma.LS.2u %*%t(R2u)
Theta2.LS.v <- R2v%*% Sigma.LS.2v %*%t(R2v)
Theta2.LS.u <- kronecker(t(RU2),diag(d2))%*%Theta2.LS.u
Theta2.LS.v <- kronecker(t(RV2),diag(d2))%*%Theta2.LS.v
return(list("Sigma"=Sigma.LS,"Theta1.u"=Theta1.LS.u,"Theta1.v"=Theta1.LS.v,
"Theta2.u"=Theta2.LS.u,"Theta2.v"=Theta2.LS.v))
}

MAR1.RRCC.SE <- function(A1,A2,k1,k2,Sigma1,Sigma2,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100){
# canonical correlation analysis
# X_t = A1 X_{t-1} A2^T + E_t
# RU1,RV1,RU2,RV2: rotation of U1,V1,U2,V2, e.g., new_U1=U1 RU1
# Cov(vec(E_t)) = Sigma.e=Sigma2 \otimes \Sigma1 : of dimension d1d2\times d1d2
# truncate vec(X_t) to mpower term, i.e. VMA(mpower)
# return Sigma.CC: asymptotic covariance matrix of (vec(\hat A1),vec(\hat A2^T))
#  Theta1.CC.u, Theta1.CC.v: asymptotic covariance matrix of vec(\hat U1), vec(\hat V1)
#  Theta2.CC.u, Theta2.CC.v: asymptotic covariance matrix of vec(\hat U2), vec(\hat V2)
d1 <- dim(A1)[1]
d2 <- dim(A2)[1]
Sigma1.inv=solve(Sigma1)
Sigma2.inv=solve(Sigma2)
Sigma.e=kronecker(Sigma2,Sigma1)
Sigma.e.inv=kronecker(Sigma2.inv,Sigma1.inv)
Sigma.ee <- array(Sigma.e,c(d1,d2,d1,d2))
B <- kronecker(A2,A1)
#Sigma.x.vec, Sigma.xx: covariance matrix of vec(X_t)
Sigma.x.vec <- Sigma.e
for(i in 1:mpower){
Sigma.x.vec <- Sigma.x.vec+(B%^%i)%*%Sigma.e%*%t(B%^%i)
}
Sigma.xx <- array(Sigma.x.vec,c(d1,d2,d1,d2))
Gamma1 <- tensor(Sigma.xx,t(A1)%*%Sigma1.inv%*%A1,c(1,3),c(1,2))
Gamma2 <- tensor(Sigma.xx,t(A2)%*%Sigma2.inv%*%A2,c(2,4),c(1,2))
Gamma3 <- kronecker(diag(1,d1),Sigma1.inv%*%A1)%*%matrix(aperm(Sigma.xx,c(3,1,4,2)),d1^2,d2^2)%*%
kronecker(t(A2)%*%Sigma2.inv,diag(1,d2))
H <- matrix(NA,d1^2+d2^2,d1^2+d2^2)
H[1:d1^2,1:d1^2] <- as.vector(A1)%*%t(as.vector(A1)) + kronecker(Gamma2,Sigma1.inv)
H[1:d1^2,(d1^2+1):(d1^2+d2^2)] <- Gamma3
H[(d1^2+1):(d1^2+d2^2),1:d1^2] <- t(Gamma3)
H[(d1^2+1):(d1^2+d2^2),(d1^2+1):(d1^2+d2^2)] <- kronecker(Sigma2.inv,Gamma1)
H.inv <- solve(H)
D1 <- Gamma2%*%t(A1)%*% ginv(A1%*%Gamma2%*%t(A1)) %*%A1
D2 <- Gamma1%*%t(A2)%*% ginv(A2%*%Gamma1%*%t(A2)) %*%A2
P1 <- Sigma1.inv%*%A1%*%ginv(t(A1)%*%Sigma1.inv%*%A1)%*%t(A1)
P2 <- Sigma2.inv%*%A2%*%ginv(t(A2)%*%Sigma2.inv%*%A2)%*%t(A2)
Sigma.Q <- matrix(NA,d1^2+d2^2,d1^2+d2^2)
M1 <- kronecker(t(A2),P1)%*%Sigma.e.inv%*%kronecker(A2,P1)
M2 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e.inv%*%kronecker(A2,P1)
M3 <- kronecker(t(A2),P1)%*%Sigma.e.inv%*%kronecker(A2,diag(d1)-P1)
M4 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e.inv%*%kronecker(A2,diag(d1)-P1)
Sigma.Q1 <- Sigma.Q2 <- Sigma.Q3 <- Sigma.Q4 <- matrix(0,d1^2,d1^2)
for(i in 1:d1){
for(j in 1:d1){
for(k in 1:d1){
for(l in 1:d1){
Sigma.Q1[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M1[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
Sigma.Q2[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M2[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
Sigma.Q3[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M3[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
Sigma.Q4[(i-1)*d1+j,(k-1)*d1+l] <- sum(Sigma.xx[i,,k,] * M4[(1:d2-1)*d1+j,(1:d2-1)*d1+l])
}
}
}
}
Sigma.Q[1:d1^2,1:d1^2] <- Sigma.Q1+kronecker(D1,diag(d1))%*%Sigma.Q2+Sigma.Q3%*%kronecker(t(D1),diag(d1))+
kronecker(D1,diag(d1))%*%Sigma.Q4%*%kronecker(t(D1),diag(d1))
M1 <- kronecker(P2,t(A1))%*%Sigma.e.inv%*%kronecker(P2,A1)
M2 <- kronecker(P2,t(A1))%*%Sigma.e.inv%*%kronecker(diag(d2)-P2,A1)
M3 <- kronecker(diag(d2)-P2,t(A1))%*%Sigma.e.inv%*%kronecker(P2,A1)
M4 <- kronecker(diag(d2)-P2,t(A1))%*%Sigma.e.inv%*%kronecker(diag(d2)-P2,A1)
Sigma.Q1 <- Sigma.Q2 <- Sigma.Q3 <- Sigma.Q4 <- matrix(0,d2^2,d2^2)
for(i in 1:d2){
for(j in 1:d2){
for(k in 1:d2){
for(l in 1:d2){
Sigma.Q1[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M1[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
Sigma.Q2[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M2[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
Sigma.Q3[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M3[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
Sigma.Q4[(i-1)*d2+j,(k-1)*d2+l] <- sum(Sigma.xx[,j,,l] * M4[(i-1)*d1+(1:d1),(k-1)*d1+(1:d1)])
}
}
}
}
Sigma.Q[(d1^2+1):(d1^2+d2^2),(d1^2+1):(d1^2+d2^2)] <- Sigma.Q1+Sigma.Q2%*%kronecker(diag(d2),t(D2))+
kronecker(diag(d2),D2)%*%Sigma.Q3+kronecker(diag(d2),D2)%*%Sigma.Q4%*%kronecker(diag(d2),t(D2))
M1 <- kronecker(t(A2),P1)%*%Sigma.e.inv%*%kronecker(P2,A1)
M2 <- kronecker(t(A2),P1)%*%Sigma.e.inv%*%kronecker(diag(d2)-P2,A1)
M3 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e.inv%*%kronecker(P2,A1)
M4 <- kronecker(t(A2),diag(d1)-P1)%*%Sigma.e.inv%*%kronecker(diag(d2)-P2,A1)
Sigma.Q1 <- Sigma.Q2 <- Sigma.Q3 <- Sigma.Q4 <- matrix(0,d1^2,d2^2)
for(i in 1:d1){
for(j in 1:d1){
for(k in 1:d2){
for(l in 1:d2){
Sigma.Q1[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M1[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
Sigma.Q2[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M2[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
Sigma.Q3[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M3[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
Sigma.Q4[(i-1)*d1+j,(k-1)*d2+l] <- sum(Sigma.xx[i,,,l] * M4[(1:d2-1)*d1+j,(k-1)*d1+(1:d1)])
}
}
}
}
Sigma.Q[1:d1^2,(d1^2+1):(d1^2+d2^2)] <- Sigma.Q1+Sigma.Q2%*%kronecker(diag(d2),t(D2))+
kronecker(D1,diag(d1))%*%Sigma.Q3+kronecker(D1,diag(d1))%*%Sigma.Q4%*%kronecker(diag(d2),t(D2))
Sigma.Q[(d1^2+1):(d1^2+d2^2),1:d1^2] <- t(Sigma.Q[1:d1^2,(d1^2+1):(d1^2+d2^2)])
Sigma.CC <- H.inv%*%Sigma.Q%*%H.inv
Sigma.CC11 <- Sigma.CC[1:d1^2,1:d1^2]
Sigma.CC11.t <- array(Sigma.CC11,c(d1,d1,d1,d1))
Sigma.CC11.t <- aperm(Sigma.CC11.t,c(2,1,4,3))
Sigma.CC11.t <- matrix(Sigma.CC11.t, d1^2, d1^2)
Sigma.CC22.t <- Sigma.CC[(d1^2+1):(d1^2+d2^2),(d1^2+1):(d1^2+d2^2)]
Sigma.CC22 <- array(Sigma.CC22.t,c(d2,d2,d2,d2))
Sigma.CC22 <- aperm(Sigma.CC22,c(2,1,4,3))
Sigma.CC22 <- matrix(Sigma.CC22, d2^2, d2^2)
svd.A1 <- svd(A1)
#k1 <- length(svd.A1$d[svd.A1$d>1e-10])
D1 <- diag(c(svd.A1$d[1:k1],1))[1:k1,1:k1] U1 <- svd.A1$u[,1:k1]
V1 <- svd.A1$v[,1:k1] if(k1<d1){ U1c <- svd.A1$u[,(k1+1):d1]
V1c <- svd.A1$v[,(k1+1):d1] }else if(k1==d1){ U1c <- 0 V1c <- 0 } e1 <- diag(d1) J1 <- matrix(0,d1^2,d1^2) for(i in 1:d1){ J1[,((i-1)*d1+1):(i*d1)] <- kronecker(diag(d1),e1[,i]) } C1 <- kronecker(A1,diag(d1)) + kronecker(diag(d1),A1)%*%J1 Sigma.CC.1u <- C1%*%Sigma.CC11%*%t(C1) C1.t <- kronecker(t(A1),diag(d1)) + kronecker(diag(d1),t(A1))%*%J1 Sigma.CC.1v <- C1.t%*%Sigma.CC11.t%*%t(C1.t) e1 <- diag(k1) L1 <- matrix(0,k1^2,k1) for(i in 1:k1){ L1[,i] <- kronecker(e1[,i],e1[,i]) } if(k1<d1){ R1u <- cbind(kronecker(diag(k1),U1), kronecker(diag(k1),U1c)) %*% rbind(solve(kronecker(D1^2,diag(k1)) - kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(U1),t(U1)), kronecker(diag(1/c(svd.A1$d[1:k1],1))[1:k1,1:k1]^2%*%t(U1), t(U1c))  )
R1v <- cbind(kronecker(diag(k1),V1), kronecker(diag(k1),V1c)) %*% rbind(solve(kronecker(D1^2,diag(k1)) -
kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(V1),t(V1)),
kronecker(diag(1/c(svd.A1$d[1:k1],1))[1:k1,1:k1]^2%*%t(V1), t(V1c)) ) }else if(k1==d1){ R1u <- kronecker(diag(k1),U1) %*% solve(kronecker(D1^2,diag(k1)) - kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(U1),t(U1)) R1v <- kronecker(diag(k1),V1) %*% solve(kronecker(D1^2,diag(k1)) - kronecker(diag(k1),D1^2)+L1%*%t(L1)) %*% (diag(k1^2)-L1%*%t(L1)) %*% kronecker(t(V1),t(V1)) } Theta1.CC.u <- R1u%*% Sigma.CC.1u %*%t(R1u) Theta1.CC.v <- R1v%*% Sigma.CC.1v %*%t(R1v) Theta1.CC.u <- kronecker(t(RU1),diag(d1))%*%Theta1.CC.u Theta1.CC.v <- kronecker(t(RV1),diag(d1))%*%Theta1.CC.v svd.A2 <- svd(A2) #k2 <- length(svd.A2$d[svd.A2$d>1e-10]) D2 <- diag(c(svd.A2$d[1:k2],1))[1:k2,1:k2]
U2 <- svd.A2$u[,1:k2] V2 <- svd.A2$v[,1:k2]
if(k1<d1){
U2c <- svd.A2$u[,(k2+1):d2] V2c <- svd.A2$v[,(k2+1):d2]
}else if(k2==d2){
U2c <- 0
V2c <- 0
}
e2 <- diag(d2)
J2 <- matrix(0,d2^2,d2^2)
for(i in 1:d2){
J2[,((i-1)*d2+1):(i*d2)] <- kronecker(diag(d2),e2[,i])
}
C2 <- kronecker(A2,diag(d2)) + kronecker(diag(d2),A2)%*%J2
Sigma.CC.2u <- C2%*%Sigma.CC22%*%t(C2)
C2.t <- kronecker(t(A2),diag(d2)) + kronecker(diag(d2),t(A2))%*%J2
Sigma.CC.2v <- C2.t%*%Sigma.CC22.t%*%t(C2.t)
e2 <- diag(k2)
L2 <- matrix(0,k2^2,k2)
for(i in 1:k2){
L2[,i] <- kronecker(e2[,i],e2[,i])
}
if(k2<d2){
R2u <- cbind(kronecker(diag(k2),U2), kronecker(diag(k2),U2c)) %*% rbind(solve(kronecker(D2^2,diag(k2)) -
kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(U2),t(U2)),
kronecker(diag(1/c(svd.A2$d[1:k2],1))[1:k2,1:k2]^2%*%t(U2), t(U2c)) ) R2v <- cbind(kronecker(diag(k2),V2), kronecker(diag(k2),V2c)) %*% rbind(solve(kronecker(D2^2,diag(k2)) - kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(V2),t(V2)), kronecker(diag(1/c(svd.A2$d[1:k2],1))[1:k2,1:k2]^2%*%t(V2), t(V2c))  )
}else if(k2==d2){
R2u <- kronecker(diag(k2),U2) %*% solve(kronecker(D2^2,diag(k2)) -
kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(U2),t(U2))
R2v <- kronecker(diag(k2),V2) %*% solve(kronecker(D2^2,diag(k2)) -
kronecker(diag(k2),D2^2)+L2%*%t(L2)) %*% (diag(k2^2)-L2%*%t(L2)) %*% kronecker(t(V2),t(V2))
}
Theta2.CC.u <- R2u%*% Sigma.CC.2u %*%t(R2u)
Theta2.CC.v <- R2v%*% Sigma.CC.2v %*%t(R2v)
Theta2.CC.u <- kronecker(t(RU2),diag(d2))%*%Theta2.CC.u
Theta2.CC.v <- kronecker(t(RV2),diag(d2))%*%Theta2.CC.v
return(list("Sigma"=Sigma.CC,"Theta1.u"=Theta1.CC.u,"Theta1.v"=Theta1.CC.v,
"Theta2.u"=Theta2.CC.u,"Theta2.v"=Theta2.CC.v))
}

MAR1.PROJ <- function(xx){
# xx: T * p * q
# X_t = LL X_{t-1} RR + E_t
# Sig = cov(vec(E_t))
# one-step projection estimation
# Return LL, RR, and estimate of Sig
dd <- dim(xx)
T <- dd[1]
p <- dd[2]
q <- dd[3]
xx.mat <- matrix(xx,T,p*q)
kroneck <- t(xx.mat[2:T,]) %*% xx.mat[1:(T-1),] %*% solve(t(xx.mat[1:(T-1),]) %*% xx.mat[1:(T-1),])
ans.projection <- projection(kroneck,r=1,q,p,q,p)
a <- svd(ans.projection$C,nu=0,nv=0)$d[1]
LL <- ans.projection$C / a RR <- t(ans.projection$B) * a
res=xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR,3,1),LL,2,2),c(1,3,2))
Sig <- matrix(tensor(res,res,1,1),p*q)/(T-1)
sd <- MAR.SE(xx, t(RR), LL, Sig)
return(list(LL=LL,RR=RR,res=res,Sig=Sig, sd=sd))
}

MAR1.LS <- function(xx,niter=50,tol=1e-6,print.true = FALSE){
# xx: T * p * q
# X_t = LL X_{t-1} RR + E_t
# Sig = cov(vec(E_t))
# LS criterion
# iterative algorithm between LL <--> RR
# Return LL, RR, and estimate of Sig
dd=dim(xx)
T <- dd[1]
p <- dd[2]
q <- dd[3]
LL.old <- diag(p)
RR.old <- diag(q)
dis <- 1
iiter <- 1

while(iiter <= niter & dis >= tol){
# estimate RR0
temp <- tensor(xx[1:(T-1),,,drop=FALSE],LL.old,2,2)  # (T-1) * q * p
AA <- tensor(temp,temp,c(1,3),c(1,3))
BB <- tensor(temp,xx[2:T,,,drop=FALSE],c(1,3),c(1,2))
RR <- solve(AA,BB)
# estimate LL0
temp <- tensor(xx[1:(T-1),,,drop=FALSE],RR,3,1)  # (T-1) * p * q
AA <- tensor(temp,temp,c(1,3),c(1,3))
BB <- t(tensor(temp,xx[2:T,,,drop=FALSE],c(1,3),c(1,3)))
LL <- t(solve(t(AA),t(BB)))
a <- svd(LL,nu=0,nv=0)$d[1] LL <- LL / a RR <- RR * a # update for the next iteration dis <- sqrt(sum((kronecker(t(RR),LL)-kronecker(t(RR.old),LL.old))^2)) LL.old <- LL RR.old <- RR iiter <- iiter + 1 if(print.true==TRUE){ print(LL) print(RR) } } res=xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR,3,1),LL,2,2),c(1,3,2)) Sig <- matrix(tensor(res,res,1,1),p*q)/(T-1) sd <- MAR.SE(xx, t(RR), LL, Sig) return(list(A1=LL,A2=t(RR),res=res,Sig=Sig,niter=iiter,sd=sd)) } MAR1.MLE <- function(xx,LL.init=NULL,Sigl.init=NULL,Sigr.init=NULL,niter=50,tol=1e-6,print.true = FALSE){ # xx: T * p * q # X_t = LL X_{t-1} RR + E_t # Sig = cov(vec(E_t)) = Sigr otimes Sigl # optimization criterion is likelihood # iterative algorithm between LL <--> RR <--> Sig_r <--> Sig_l # Return LL, RR, Sigl, Sigr dd=dim(xx) T <- dd[1] p <- dd[2] q <- dd[3] if(is.null(LL.init)){ LL.old <- diag(p) } else{ LL.old <- LL.init } RR.old <- diag(q) if(is.null(Sigl.init)){ Sigl.old <- diag(p) } else{ Sigl.old <- Sigl.init } if(is.null(Sigr.init)){ Sigr.old <- diag(q) } else{ Sigr.old <- Sigr.init } dis <- 1 iiter <- 1 while(iiter <= niter & dis >= tol){ Sigl.inv.old <- ginv(Sigl.old) Sigr.inv.old <- ginv(Sigr.old) # estimate RR0 temp1 <- tensor(xx[1:(T-1),,,drop=FALSE],LL.old,2,2) # (T-1) * q * p temp2 <- tensor(temp1,Sigl.inv.old,3,2) # (T-1) * q * p AA <- tensor(temp1,temp2,c(1,3),c(1,3)) BB <- tensor(temp2,xx[2:T,,,drop=FALSE],c(1,3),c(1,2)) RR <- solve(AA,BB) # estimate LL0 temp1 <- tensor(xx[1:(T-1),,,drop=FALSE],RR,3,1) # (T-1) * p * q temp2 <- tensor(temp1,Sigr.inv.old,3,1) # (T-1) * p * q AA <- tensor(temp1,temp2,c(1,3),c(1,3)) BB <- t(tensor(temp2,xx[2:T,,,drop=FALSE],c(1,3),c(1,3))) LL <- t(solve(t(AA),t(BB))) res=xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR,3,1),LL,2,2),c(1,3,2)) temp <- tensor(res,Sigl.inv.old,2,1) # (T-1) * q * p Sigr <- tensor(temp,res,c(1,3),c(1,2))/T/p temp <- tensor(res,ginv(Sigr),3,1) # (T-1) * p * q Sigl <- tensor(temp,res,c(1,3),c(1,3))/T/q a <- svd(LL,nu=0,nv=0)$d[1]
LL <- LL / a
RR <- RR * a
a <- eigen(Sigl)$values[1] Sigl <- Sigl / a Sigr <- Sigr * a # update for the next iteration dis1 <- sqrt(sum((kronecker(t(RR),LL)-kronecker(t(RR.old),LL.old))^2)) dis2 <- sqrt(sum((kronecker(Sigr,Sigl)-kronecker(Sigr.old,Sigl.old))^2)) dis <- max(dis1,dis2) LL.old <- LL RR.old <- RR Sigr.old <- Sigr Sigl.old <- Sigl iiter <- iiter + 1 if(print.true==TRUE){ print(LL) print(RR) } } Sig <- kronecker(Sigr,Sigl) sd <- MAR.SE(xx, t(RR), LL, Sig) return(list(A1=LL,A2=t(RR),res=res,Sig1=Sigl,Sig2=Sigr,Sig=Sig,niter=iiter,sd=sd)) } MAR.SE <- function(xx, B, A, Sigma){ dd <- dim(xx) T <- dd[1] p <- dd[2] q <- dd[3] BX <- tensor(xx,B,3,2) # Tpq AX <- tensor(xx,A,2,2) # Tqp BXI <- apply(BX,1,function(x){kronecker(t(x),diag(p))}) AXI <- apply(AX,1,function(x){kronecker(diag(q),t(x))}) BXI.array <- array(BXI,c(p*q,p^2,T)) #pq*p^2*T AXI.array <- array(AXI,c(p*q,q^2,T)) #pq*q^2*T # W transpose Wt <- abind(BXI.array,AXI.array,along=2) #pq*(p^2+q^2)*T EWWt <- tensor(Wt,Wt,c(1,3),c(1,3))/T #(p^2+q^2)*(p^2+q^2) alpha <- as.vector(A) beta <- as.vector(t(B)) gamma <- c(alpha,rep(0,q^2)) H <- EWWt + gamma %*% t(gamma) #(p^2+q^2)*(p^2+q^2) WSigma <- tensor(Wt,Sigma,1,1) #(p^2+q^2)*T*pq EWSigmaWt <- tensor(WSigma,Wt,c(3,2),c(1,3))/T Hinv <- solve(H) Xi <- Hinv %*% EWSigmaWt %*% Hinv Xi <- Xi/T } tenAR.VAR <- function(xx, P){ dd=dim(xx) t <- dd[1] n <- prod(dd[-1]) if (prod(dd[-1]) > t){stop("sample size T too small")} yy=apply(xx, MARGIN=1, as.vector) x = 0 for (l in c(1:P)){ if (l == 1){x = yy[,(P):(t-1)]} else {x = rbind(x, yy[,(P+1-l):(t-l)])} } y = yy[,(P+1):t] coef = t(solve(x %*% t(x)) %*% x %*% t(y)) res= y - coef %*% x phi = list() for (l in c(1:P)){ phi[[l]] = coef[,(n*(l-1)+1):(n*l)] } return(list(coef=phi, res=res)) } tenAR.A <- function(dim,R,P,rho){ K <- length(dim) A <- lapply(1:P, function(p) {lapply(1:max(1,R[p]), function(j) {lapply(1:K, function(i) {diag(dim[i])})})}) for (p in c(1:P)){ if (R[p] == 0) next for (j in c(1:R[p])){ for (i in c(1:K)){ A[[p]][[j]][[i]] <- matrix(rnorm(dim[i]^2), c(dim[i],dim[i])) } } } eigen = M.eigen(A, R, P, dim) for (l in c(1:P)){ if (R[l] == 0) next for (j in c(1:R[l])){ A[[l]][[j]][[K]] <- rho * A[[l]][[j]][[K]]/ eigen } A[[l]] <- fro.order(fro.rescale(A[[l]])) } eigen = M.eigen(A, R, P, dim) return(A) } tenAR.PROJ <- function(xx,R,P){ dim <- dim(xx)[-1] mm <- tenAR.VAR(xx, P)$coef
A = list()
for (p in c(1:P)){
if (is.na(R[p])) stop("p != length(R)")
if (R[p] == 0) next
tt <- aperm(trearrange(mm[[p]], rev(dim)))
A[[p]] <- fro.order(fro.rescale(ten.proj(tt, dim, R[p])))
}
return(list(A=A))
}

tenAR.LS <- function(xx, R, P, init.A=NULL, niter=150, tol=1e-5,print.true = FALSE){
if (!(mode(xx) == "S4")) {xx <- as.tensor(xx)}
dim <- dim(xx)[-1]
K <- length(dim)
t <- dim(xx)[1]
if (K==2){

if (is.null(init.A)) {

A.old = list()
for (p in c(1:P)){
if (is.na(R[p])) stop("p != length(R)")
if (R[p] == 0) next
A.old[[p]] <- lapply(1:R[p], function(j) {lapply(1:K, function(i) {0.5*diag(dim[i])})})
}

} else {A.old <- init.A}
}
if (K==3) {if (is.null(init.A)) {A.old <- tenAR.PROJ(xx@data,R,P)$A} else {A.old <- init.A}} A.new <- A.old Tol <- tol*sqrt(sum(dim^2))*sum(R) dis <- 1 iiter <- 1 while(iiter <= niter & dis >= tol){ dis3 <- 0 for (p in c(1:P)){ if (R[p] == 0) next # print(p) for (r in c(1:R[p])){ # print(r) for (k in c(K:1)){ # update last matrix first # print(k) # tic("step 1") # temp <- tl(xx, A.new[[p]][[r]][-k], k)[(1+P-p):(t-p),,,,drop=FALSE] temp <- rTensor::ttl(xx, A.new[[p]][[r]][-k], c(2:(K+1))[-k]) temp <- myslice(temp,K,1+P-p,t-p) L1 <- 0 # toc() # tic("step 2") for (l in c(1:P)){ if (R[l] == 0) next if (l == p){if (R[l] > 1){L1 <- L1 + Reduce("+",lapply(c(1:R[l])[-r], function(n) {rTensor::ttl(myslice(xx, K, 1+P-l, t-l), A.new[[l]][[n]], (c(1:K) + 1))}))} } else {L1 <- L1 + Reduce("+",lapply(c(1:R[l]), function(n) {rTensor::ttl(myslice(xx, K, 1+P-l, t-l), A.new[[l]][[n]], (c(1:K) + 1))}))} } temp2 <- myslice(xx, K, 1+P, t) - L1 # toc() # tic("step 3") RR <- tensor(temp@data,temp@data,c(1:(K+1))[-(k+1)],c(1:(K+1))[-(k+1)]) # toc() # tic("step 4") LL <- tensor(temp2@data,temp@data,c(1:(K+1))[-(k+1)],c(1:(K+1))[-(k+1)]) # toc() # tic("step 5") A.new[[p]][[r]][[k]] <- LL %*% ginv(RR) # toc() dis3 <- dis3 + min(sum((A.new[[p]][[r]][[k]] - A.old[[p]][[r]][[k]])^2), sum((-A.new[[p]][[r]][[k]] - A.old[[p]][[r]][[k]])^2)) } } } # print("done once") # tic("svd") for (p in c(1:P)){ if (R[p] == 0) next A.new[[p]] <- svd.rescale(A.new[[p]]) } # toc() dis <- sqrt(dis3) A.old <- A.new iiter <- iiter + 1 if (print.true == TRUE){ print(dis) print(paste('iiter num=',iiter)) } } for (p in c(1:P)){ if (R[p] == 0) next A.new[[p]] <- fro.order(fro.rescale(A.new[[p]])) } res <- ten.res(xx,A.new,P,R,K,t)@data Sig <- matrix(tensor(res,res,1,1),prod(dim))/(t-1) if (K==3 & P==1){ cov = tenAR.SE.LSE(xx, A.new[[1]], Sig) sd <- covtosd(cov, dim, R)} # temporarily for P=1 only} bic <- IC(xx, res, R, t, dim) return(list(A=A.new,niter=iiter,Sig=Sig,res=res,cov=cov,sd=sd,BIC=bic)) } tenAR.MLE <- function(xx, R, P, init.A=NULL, init.sig=NULL, niter=150,tol=1e-5, print.true = FALSE){ if (!(mode(xx) == "S4")) {xx <- as.tensor(xx)} dim <- dim(xx)[-1] K <- length(dim) t <- dim(xx)[1] if (is.null(init.sig)) {Sig.old <- lapply(1:K, function(i) {diag(dim[i])})} else {Sig.old <- init.sig} Sig.new <- Sig.old Sig.new.inv <- lapply(1:K, function (k) {solve(Sig.new[[k]])}) if (K >= 3){if (is.null(init.A)) {A.old <- tenAR.PROJ(xx@data,R,P)$A} else {A.old <- init.A}}
if (K==2){

if (is.null(init.A)) {

A.old = list()
for (p in c(1:P)){
if (is.na(R[p])) stop("p != length(R)")
if (R[p] == 0) next
A.old[[p]] <- lapply(1:R[p], function(j) {lapply(1:K, function(i) {0.5*diag(dim[i])})})
}

} else {A.old <- init.A}
}
A.new <- A.old
Tol <- tol*sqrt(sum(dim^2))*P*sum(R)
dis <- 1
iiter <- 1
while(iiter <= niter & dis >= Tol){
stopifnot(dis <= 1e3)
dis3 <- 0
for (p in c(1:P)){
for (r in c(1:R[p])){
for (k in c(K:1)){
res.old <- ten.res(xx,A.new,P,R,K,t)
rs <- rTensor::ttl(res.old, Sig.new.inv[-k], c(2:(K+1))[-k])
Sig.new[[k]] <- tensor(res.old@data, rs@data, c(1:(K+1))[-(k+1)],c(1:(K+1))[-(k+1)])/(t-1)/prod(dim[-k])
Sig.new.inv <- lapply(1:K, function (k) {ginv(Sig.new[[k]])})
}
for (k in c(K:1)){
sphi <-  lapply(1:K, function (k) {Sig.new.inv[[k]] %*% (A.new[[p]][[r]][[k]])})
temp <- ttl(xx, A.new[[p]][[r]][-k], c(2:(K+1))[-k])
temp = myslice(temp, K, 1+P-p, t-p)
temp1 <- ttl(xx, sphi[-k],c(2:(K+1))[-k])
temp1 = myslice(temp1, K, 1+P-p, t-p)
L1 <- 0
for (l in c(1:P)){
if (l == p){
if (R[l] > 1){
L1 <- L1 + Reduce("+",lapply(c(1:R[l])[-r], function(n) {rTensor::ttl(myslice(xx, K, 1+P-l, t-l), A.new[[l]][[n]], (c(1:K) + 1))}))
}
} else {
L1 <- L1 + Reduce("+",lapply(c(1:R[l]), function(n) {rTensor::ttl(myslice(xx, K, 1+P-l, t-l), A.new[[l]][[n]], (c(1:K) + 1))}))
}
}
temp2 <-  myslice(xx,K,1+P,t) - L1
RR <- tensor(temp@data,temp1@data,c(1:(K+1))[-(k+1)],c(1:(K+1))[-(k+1)])
LL <- tensor(temp2@data,temp1@data,c(1:(K+1))[-(k+1)],c(1:(K+1))[-(k+1)])
A.new[[p]][[r]][[k]] <- LL %*% ginv(RR)
dis3 <- dis3 + min(sum((A.new[[p]][[r]][[k]] - A.old[[p]][[r]][[k]])^2), sum((-A.new[[p]][[r]][[k]] - A.old[[p]][[r]][[k]])^2))
}
}
A.new[[p]] <- svd.rescale(A.new[[p]])
Sig.new <- eigen.rescale(list(Sig.new))[[1]]
}
dis <- sqrt(dis3)
Sig.old <- Sig.new
A.old <- A.new
iiter <- iiter + 1
if (print.true == TRUE){
print(dis)
print(paste('iiter num=',iiter))
}
}
for (p in c(1:P)){
A.new[[p]] <- fro.order(fro.rescale(A.new[[p]]))
}
res <- ten.res(xx,A.new,P,R,K,t)@data
Sig <- matrix(tensor(res,res,1,1),prod(dim))/(t-1)
if (K==3 & P==1){
cov <- tenAR.SE.MLE(xx@data, A.new[[1]], Sig) # temporarily for P=1 only
sd <- covtosd(cov, dim, R)
}
bic <- IC(xx, res, R, t, dim)
return(list(A=A.new, SIGMA=Sig.new, niter=iiter, Sig=Sig, res=res, cov=cov, sd=sd, BIC=bic))
}

MAR1.RR <- function(xx, k1, k2, niter=200, tol=1e-4, A1.init=NULL, A2.init=NULL, print.true=FALSE){
# xx: T * p * q
# X_t = LL X_{t-1} RR' + E_t ### NOTE! This function is written with RR'.
# Sig = cov(vec(E_t)) = Sigr \otimes Sigl
# optimization criterion is likelihood
# iterative algorithm between LL <--> Sigma_l <--> RR <--> Sig_r
# Return LL, RR, Sigl, Sigr
dd=dim(xx)
T <- dd[1]
p <- dd[2]
q <- dd[3]
if(is.null(A1.init)){
LL.old <- diag(p)
} else{
LL.old <- A1.init
}
if(is.null(A2.init)){
RR.old <- diag(q)
} else{
RR.old <- A2.init
}
Tol=tol*sqrt(p^2+q^2)
dis <- 1
iiter <- 1

while(iiter <= niter & dis >= Tol){
## Save old
LL.oold=LL.old
RR.oold=RR.old

# estimate LL0
temp1 <- tensor(xx[1:(T-1),,,drop=FALSE],RR.old,3,2)  # (T-1) * p * q
AA <- tensor(temp1,temp1,c(1,3),c(1,3))
BB <- tensor(xx[2:T,,,drop=FALSE],temp1,c(1,3),c(1,3))
LL <- BB%*%ginv(AA)
U <- svd(LL%*%t(BB))$u[,1:k1] LL <- U%*%t(U)%*%LL # update for next iteration a=svd(LL,nu=0,nv=0)$d[1]
LL=LL/a
dis3=sum((LL-LL.old)^2)
LL.old <- LL

# estimate RR0
temp1 <- tensor(xx[1:(T-1),,,drop=FALSE],LL.old,2,2)  # (T-1) * q * p
AA <- tensor(temp1,temp1,c(1,3),c(1,3))
BB <- tensor(xx[2:T,,,drop=FALSE],temp1,c(1,2),c(1,3))
RR <- BB%*%ginv(AA)
U <- svd(RR%*%t(BB))$u[,1:k2] RR <- U%*%t(U)%*%RR # update for next iteration dis3=dis3+sum((RR-RR.old)^2) RR.old <- RR # update for the next iteration dis1 <- sqrt(sum((kronecker(t(RR),LL)-kronecker(t(RR.oold),LL.oold))^2)) dis3 = sqrt(dis3) dis <- dis3 iiter <- iiter + 1 #print(max(abs(eigen(LL)$values)))
#print(max(abs(eigen(RR)$values))) if(print.true==TRUE){ print(dis) print(paste('iiter num=',iiter)) } } a <- sqrt(sum(LL^2)) LL <- LL / a RR <- RR * a res=xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR,3,2),LL,2,2),c(1,3,2)) Sig <- matrix(tensor(res,res,1,1),p*q)/(T-1) bic <- T*p*q*log(sum(res^2/(T*p*q))) ##+log(T*p*q)*(bic.penalty(p,k1)+bic.penalty(q,k2)) cov <- matAR.RR.se(LL,RR,k1,k2,method="RRLSE",Sigma.e=Sig,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100)$Sigma
sd <- list(array(diag(cov)[1:p^2], c(p,p)), as.array(diag(cov)[(p^2+1):(p^2+q^2)], c(q,q)))
return(list(A1=LL,A2=RR,res=res,Sig=Sig,BIC=bic,niter=iiter-1,cov=cov,sd=sd))
}

MAR1.CC <- function(xx,k1,k2,A1.init=NULL,A2.init=NULL,Sigl.init=NULL,Sigr.init=NULL,niter=200,tol=1e-4,print.true = FALSE){
# xx: T * p * q
# X_t = LL X_{t-1} RR' + E_t ### NOTE! This function is written with RR'.
# Sig = cov(vec(E_t)) = Sigr \otimes Sigl
# optimization criterion is likelihood
# iterative algorithm between LL <--> Sigma_l <--> RR <--> Sig_r
# Return LL, RR, Sigl, Sigr
dd=dim(xx)
T <- dd[1]
p <- dd[2]
q <- dd[3]
if(is.null(A1.init)){
LL.old <- diag(p)
} else{
LL.old <- A1.init
}
if(is.null(A2.init)){
RR.old <- diag(q)
} else{
RR.old <- A2.init
}
if(is.null(Sigl.init)){
Sigl.old <- diag(p)
} else{
Sigl.old <- Sigl.init
}
if(is.null(Sigr.init)){
Sigr.old <- diag(q)
} else{
Sigr.old <- Sigr.init
}
Tol=tol*sqrt(p^2+q^2)
dis <- 1
iiter <- 1

while(iiter <= niter & dis >= Tol){

## Save old
LL.oold=LL.old
RR.oold=RR.old
Sigl.oold=Sigl.old
Sigr.oold=Sigr.old

# estimate LL0 and Sigl
Sigr.inv.old <- ginv(Sigr.old)
temp1 <- tensor(xx[1:(T-1),,,drop=FALSE],RR.old,3,2)  # (T-1) * p * q
temp2 <- tensor(temp1,Sigr.inv.old,3,1)  # (T-1) * p * q
AA <- tensor(temp1,temp2,c(1,3),c(1,3))
BB <- tensor(xx[2:T,,,drop=FALSE],temp2,c(1,3),c(1,3))
LL <- BB%*%ginv(AA)
res <- xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR.old,3,2),LL,2,2),c(1,3,2)) # (T-1) * p * q
temp <- tensor(res,Sigr.inv.old,3,1) # (T-1) * p * q
Sigl <- tensor(temp,res,c(1,3),c(1,3))/T/q
Sigl.spec <- eigen(Sigl)
Sigl.root <- Sigl.spec$vectors%*%diag(sqrt(Sigl.spec$values))%*%t(Sigl.spec$vectors) Sigl.root.inv <- Sigl.spec$vectors%*%diag(1/sqrt(Sigl.spec$values))%*%t(Sigl.spec$vectors)
U <- svd(Sigl.root.inv%*%LL%*%t(BB)%*%Sigl.root.inv)$u[,1:k1] LL <- Sigl.root%*%U%*%t(U)%*%Sigl.root.inv%*%LL res <- xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR.old,3,2),LL,2,2),c(1,3,2)) # (T-1) * p * q temp <- tensor(res,Sigr.inv.old,3,1) # (T-1) * p * q Sigl <- tensor(temp,res,c(1,3),c(1,3))/T/q # update for next iteration a=svd(LL,nu=0,nv=0)$d[1]
LL=LL/a
dis3=sum((LL-LL.old)^2)
LL.old <- LL
Sigl.old <- Sigl

# estimate RR0 and Sigr
Sigl.inv.old <- ginv(Sigl.old)
temp1 <- tensor(xx[1:(T-1),,,drop=FALSE],LL.old,2,2)  # (T-1) * q * p
temp2 <- tensor(temp1,Sigl.inv.old,3,1)  # (T-1) * q * p
AA <- tensor(temp1,temp2,c(1,3),c(1,3))
BB <- tensor(xx[2:T,,,drop=FALSE],temp2,c(1,2),c(1,3))
RR <- BB%*%ginv(AA)
res <- xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR,3,2),LL.old,2,2),c(1,3,2)) # (T-1) * p * q
temp <- tensor(res,Sigl.inv.old,2,1) # (T-1) * q * p
Sigr <- tensor(temp,res,c(1,3),c(1,2))/T/p
Sigr.spec <- eigen(Sigr)
Sigr.root <- Sigr.spec$vectors%*%diag(sqrt(Sigr.spec$values))%*%t(Sigr.spec$vectors) Sigr.root.inv <- Sigr.spec$vectors%*%diag(1/sqrt(Sigr.spec$values))%*%t(Sigr.spec$vectors)
U <- svd(Sigr.root.inv%*%RR%*%t(BB)%*%Sigr.root.inv)$u[,1:k2] RR <- Sigr.root%*%U%*%t(U)%*%Sigr.root.inv%*%RR res <- xx[2:T,,,drop=FALSE] - aperm(tensor(tensor(xx[1:(T-1),,,drop=FALSE],RR,3,2),LL.old,2,2),c(1,3,2)) # (T-1) * p * q temp <- tensor(res,Sigl.inv.old,2,1) # (T-1) * q * p Sigr <- tensor(temp,res,c(1,3),c(1,2))/T/p # update for next iteration dis3=dis3+sum((RR-RR.old)^2) RR.old <- RR Sigr.old <- Sigr a <- eigen(Sigl)$values[1]
Sigl <- Sigl / a
Sigr <- Sigr * a
### cat(eigen(Sigl)$values[1]," ",eigen(Sigr)$values[1], " ")

# update for the next iteration
dis1 <- sqrt(sum((kronecker(t(RR),LL)-kronecker(t(RR.oold),LL.oold))^2))
### cat(dis1," ")
dis2 <- sqrt(sum((kronecker(Sigr,Sigl)-kronecker(Sigr.oold,Sigl.oold))^2))
### cat(dis2," ",dis3,"\n")
dis3 = sqrt(dis3)
### dis <- max(dis1,dis2)
dis <- dis3
Sigr.old <- Sigr
Sigl.old <- Sigl
iiter <- iiter + 1
if(print.true==TRUE){
print(LL)
print(RR)
}
}
a <- sqrt(sum(LL^2))
LL <- LL / a
RR <- RR * a
Sig <- kronecker(Sigr,Sigl)
cov <- matAR.RR.se(LL,RR,k1,k2,method="RRMLE",Sigma1=Sigl,Sigma2=Sigr,RU1=diag(k1),RV1=diag(k1),RU2=diag(k2),RV2=diag(k2),mpower=100)$Sigma sd <- list(array(diag(cov)[1:p^2], c(p,p)), as.array(diag(cov)[(p^2+1):(p^2+q^2)], c(q,q))) return(list(A1=LL,A2=RR,res=res,Sig1=Sigl,Sig2=Sigr,Sig=Sig,niter=iiter-1, cov=cov, sd=sd)) } tenAR.SE.LSE <- function(xx, A.true, Sigma){ if (mode(xx) != "S4") {xx <- rTensor::as.tensor(xx)} r <- length(A.true) dim <- xx@modes[-1] m1 <- dim[1] m2 <- dim[2] m3 <- dim[3] k <- length(dim) T <- xx@modes[[1]] ndim <- m1^2+m2^2+m3^2 Gamma <- matrix(0,r*ndim,r*ndim) r1 <- matrix(0, r*ndim, 1) for (x in c(1:r)){ for (y in c(1:(k-1))) { if (y == 1) { a <- as.matrix(as.vector(A.true[[x]][[y]])) r1[((x-1)*ndim + 1):((x-1)*ndim + m1^2),] <- a } else if (y == 2) { a <- as.matrix(as.vector(A.true[[x]][[y]])) r1[((x-1)*ndim + m1^2 + 1):((x-1)*ndim + m1^2 + m2^2),] <- a } else { stop("Only Support k=3 but now k>3") } Gamma <- Gamma + r1 %*% t(r1) } } Hdim <- r*ndim HT <- array(1,c(T,Hdim,Hdim)) WT <- array(1,c(T,Hdim,prod(dim))) #T*r(m1^2+m2^2+m^3)*(m1m2m3) for (i in c(1:r)) { C <- A.true[[i]][[3]] B <- A.true[[i]][[2]] A <- A.true[[i]][[1]] for (t in c(1:T)){ w1 <- k_unfold(xx[t,,,], m = 1)@data %*% t(kronecker(C,B)) w2 <- k_unfold(xx[t,,,], m = 2)@data %*% t(kronecker(C,A)) w3 <- k_unfold(xx[t,,,], m = 3)@data %*% t(kronecker(B,A)) w <- rbind(kronecker(w1,diag(m1)) ,kronecker(w2,diag(m2)) %*% kronecker(diag(m3),pm(m2,m1)), kronecker(w3,diag(m3)) %*% pm(m3,m2*m1)) WT[t, ((i-1)*ndim + 1):(i*ndim),] <- w } } WSigma <- tensor(WT,Sigma,3,1) #T*(m1^2+m2^2+m^3)*(m1m2m3) EWSigmaWt <- tensor(WSigma,WT,c(3,1),c(3,1))/T H <- tensor(WT,WT,c(3,1),c(3,1))/T + Gamma #r(m1^2+m2^2+m^3)*r(m1^2+m2^2+m^3) Hinv <- solve(H) Xi <- Hinv %*% EWSigmaWt %*% Hinv return(Xi) } tenAR.SE.MLE <- function(xx, A.true, Sigma){ if (mode(xx) != "S4") {xx <- rTensor::as.tensor(xx)} r <- length(A.true) dim <- xx@modes[-1] m1 <- dim[1] m2 <- dim[2] m3 <- dim[3] k <- length(dim) T <- xx@modes[[1]] ndim <- m1^2+m2^2+m3^2 Gamma <- matrix(0,r*ndim,r*ndim) r1 <- matrix(0, r*ndim, 1) for (x in c(1:r)){ for (y in c(1:(k-1))) { if (y == 1) { a <- as.matrix(as.vector(A.true[[x]][[y]])) r1[((x-1)*ndim + 1):((x-1)*ndim + m1^2),] <- a } else if (y == 2) { a <- as.matrix(as.vector(A.true[[x]][[y]])) r1[((x-1)*ndim + m1^2 + 1):((x-1)*ndim + m1^2 + m2^2),] <- a } else { stop("Only Support k=3 but now k>3") } Gamma <- Gamma + r1 %*% t(r1) } } Hdim <- r*ndim HT <- array(1,c(T,Hdim,Hdim)) WT <- array(1,c(T,Hdim,prod(dim))) #T*r(m1^2+m2^2+m^3)*(m1m2m3) for (i in c(1:r)) { C <- A.true[[i]][[3]] B <- A.true[[i]][[2]] A <- A.true[[i]][[1]] for (t in c(1:T)){ w1 <- k_unfold(xx[t,,,], m = 1)@data %*% t(kronecker(C,B)) w2 <- k_unfold(xx[t,,,], m = 2)@data %*% t(kronecker(C,A)) w3 <- k_unfold(xx[t,,,], m = 3)@data %*% t(kronecker(B,A)) w <- rbind(kronecker(w1,diag(m1)) ,kronecker(w2,diag(m2)) %*% kronecker(diag(m3),pm(m2,m1)), kronecker(w3,diag(m3)) %*% pm(m3,m2*m1)) WT[t, ((i-1)*ndim + 1):(i*ndim),] <- w } } WSigma <- tensor(WT,solve(Sigma),3,1) #T*(m1^2+m2^2+m^3)*(m1m2m3) EWSigmaWt <- tensor(WSigma,WT,c(3,1),c(3,1))/T H <- EWSigmaWt + Gamma Hinv <- solve(H) Xi <- Hinv %*% EWSigmaWt %*% Hinv return(Xi) } tenAR.bic <- function(xx, rmax=5){ if (mode(xx) != "S4") {xx <- rTensor::as.tensor(xx)} dim <- xx@modes[-1] t <- xx@modes[[1]] ans <- c() for (r in c(1:rmax)){ est <- tenAR.LS(xx, R=r, P=1) res <- est$res
ans[r] <- IC(xx,res,r,t, dim)
}
which.min(ans)
}

#' Plot Matrix-Valued Time Series
#'
#' Plot matrix-valued time series, can be also used to plot a given slice of a tensor-valued time series.
#'@name mplot
#'@rdname mplot
#'@aliases mplot
#'@export
#'@importFrom graphics par
#'@param xx  \eqn{T \times d_1 \times d_2} matrix-valued time series. Note that the number of mode is 3, where the first mode is time.
#'@return a figure.
#'@examples
#' dim <- c(3,3,3)
#' xx <- tenAR.sim(t=50, dim, R=2, P=1, rho=0.5, cov='iid')
#' mplot(xx[1:30,,,1])
mplot <- function(xx){
if (mode(xx) == "S4"){xx = xx@data}
dim = dim(xx)
time = array(c(1:dim[1]),dim[1])
opar <- par(mfrow=c(dim[2],dim[3]),mai=0.05*c(1,1,1,1),oma=c(2,2,0,0))
on.exit(par(opar))
for(i in 1:dim[2]){
for(j in 1:dim[3]){
if(i!=dim[2] & j!=1){
plot(time,xx[,i,j],type='l',xaxt='n',yaxt='n',ylim=range(xx[,i,]))
}
if(i!=dim[2] & j==1){
plot(time,xx[,i,j],type='l',xaxt='n',ylim=range(xx[,i,]))
}
if(i==dim[2] & j!=1){
plot(time,xx[,i,j],type='l',yaxt='n',ylim=range(xx[,i,]))
}
if(i==dim[2] & j==1){
plot(time,xx[,i,j],type='l',ylim=range(xx[,i,]))
}
}
}

}

#' Plot ACF of Matrix-Valued Time Series
#'
#' Plot ACF of matrix-valued time series, can be also used to plot ACF of a given slice of a tensor-valued time series.
#'@name mplot.acf
#'@rdname mplot.acf
#'@aliases mplot.acf
#'@export
#'@importFrom graphics par
#'@importFrom graphics plot
#'@importFrom stats acf
#'@param xx  \eqn{T \times d_1 \times d_2} matrix-valued time series. Note that the number of mode is 3, where the first mode is time.
#'@return a figure.
#'@examples
#' dim <- c(3,3,3)
#' xx <- tenAR.sim(t=50, dim, R=2, P=1, rho=0.5, cov='iid')
#' mplot.acf(xx[1:30,,,1])
mplot.acf <- function(xx){
if (mode(xx) == "S4"){xx = xx@data}
dim = dim(xx)
opar <- par(mfrow=c(dim[2],dim[3]),mai=0.05*c(1,1,1,1),oma=c(2,2,0,0))
on.exit(par(opar))
for(i in 1:dim[2]){
for(j in 1:dim[3]){
if(i!=dim[2] & j!=1){
acf(xx[,i,j])
}
if(i!=dim[2] & j==1){
acf(xx[,i,j])
}
if(i==dim[2] & j!=1){
acf(xx[,i,j])
}
if(i==dim[2] & j==1){
acf(xx[,i,j])
}
}
}
}

#' Predictions for Tensor Autoregressive Models
#'
#' Prediction based on the tensor autoregressive model or reduced rank MAR(1) model. If \code{rolling = TRUE}, returns the rolling forecasts.
#'@name tenAR.predict
#'@rdname tenAR.predict
#'@aliases tenAR.predict
#'@export
#'@importFrom abind abind
#'@param object a model object returned by \code{tenAR.est()}.
#'@param xx \eqn{T^{\prime} \times d_1 \times \cdots \times d_K} tensor time series.
#'@param rolling TRUE or FALSE, rolling forecast, is FALSE by default.
#'@param n0 only if \code{rolling = TRUE}, the starting point of rolling forecast.
#'@return
#'a tensor time series of length \code{n.head} if \code{rolling = FALSE};
#'
#'a tensor time series of length \eqn{T^{\prime} - n_0 - n.head + 1} if \code{rolling = TRUE}.
#'@seealso 'predict.ar' or 'predict.arima'
#'@examples
#' set.seed(333)
#' dim <- c(2,2,2)
#' t = 20
#' xx <- tenAR.sim(t, dim, R=2, P=1, rho=0.5, cov='iid')
#' est <- tenAR.est(xx, R=1, P=1, method="LSE")
#' pred <- tenAR.predict(est, xx, n.head = 1)
#' # rolling forcast
#' n0 = t - min(50,t/2)
#' pred.rolling <- tenAR.predict(est, xx, n.head = 5, rolling=TRUE, n0)
#'
#' # prediction for reduced rank MAR(1) model
#' dim <- c(2,2)
#' t = 20
#' xx <- tenAR.sim(t, dim, R=1, P=1, rho=0.5, cov='iid')
#' est <- matAR.RR.est(xx, method="RRLSE", k1=1, k2=1)
#' pred <- tenAR.predict(est, xx, n.head = 1)
#' # rolling forcast
#' n0 = t - min(50,t/2)
#' pred.rolling <- tenAR.predict(est, xx, n.head = 5, rolling=TRUE, n0)
tenAR.predict <- function(object, xx, n.head, rolling=FALSE, n0=NULL){
if (is.null(object$SIGMA)){method = "LSE"} else {method = "MLE"} if (!is.null(object$A1)){
A <- list(list(list(object$A1, object$A2)))
if (is.null(object$Sig1)){ method = "RRLSE" } else { method = "RRMLE" } } else { A <- object$A
}

if (mode(xx) != "S4") {xx <- rTensor::as.tensor(xx)}

if (rolling == TRUE){
}

P <- length(A)
R <- length(A[[1]])
K <- xx@num_modes - 1
dim <- xx@modes
for(tt in ttt){
tti <- tt - ttt[1] + 1
L1 = 0
for (l in c(1:P)){
x0 <- as.tensor(array((k_unfold(xx, m=1))@data[tt-l,], c(1,dim[-1])))
L1 <- L1 + Reduce("+",lapply(c(1:R), function(n) {(rTensor::ttl(x0, A[[l]][[n]], (c(1:K) + 1)))}))
}
xx.pred[tti, ] <- k_unfold(L1, m=1)@data
xx <- as.tensor(abind(xx@data, L1@data, along=1))
}
}

predict.rolling <- function(A, xx, n.head, method, n0){
if ((method == "RRLSE") || (method == "RRMLE")){
k1 <- rankMatrix(A[[1]][[1]][[1]])
k2 <- rankMatrix(A[[1]][[1]][[2]])
}

P <- length(A)
R <- length(A[[1]])
K <- xx@num_modes - 1
dim <- xx@modes
t <- dim[1]
if(is.null(n0)){n0 = t - min(50,t/2)}
for(tt in ttt){
tti <- tt - ttt[1] + 1
xt <- as.tensor(array((k_unfold(xx, m=1))@data[1:(tt-1),], c(tt-1,dim[-1])))
if (tti > 1){
if (method == "RRLSE"){
model <- MAR1.RR(xt@data, k1, k2)
A <- list(list(list(model$A1, model$A2)))
} else if (method == "RRMLE"){
model <- MAR1.CC(xt@data, k1, k2)
A <- list(list(list(model$A1, model$A2)))
} else {
model = tenAR.est(xt, R, P, method)
A <- model\$A
}
}
L1 = 0
for (l in c(1:P)){
x0 <- as.tensor(array((k_unfold(xt, m=1))@data[tt-l,], c(1,dim[-1])))
L1 <- L1 + Reduce("+",lapply(c(1:R), function(n) {(rTensor::ttl(x0, A[[l]][[n]], (c(1:K) + 1)))}))
}
xx.pred[tti, ] <- k_unfold(L1, m=1)@data
#xt <- as.tensor(abind(xt@data, L1@data, along=1))
}
}

.getpos <- function(mode, rank){
pos = 0
for (k in c(1:length(mode))){
if (k > 1){mode[k] = mode[k] - 1}
pos = pos + rank[k]*mode[k]
}
return(pos)
}

.getrank <- function(dim){
rank = array(1, length(dim))
for (k in c(1:length(dim))){
if (k > 1){ for (q in c(1:(k-1))){rank[k] = rank[k]*(rev(dim)[q])}}
}
return(rank)
}

.remove.mean <- function(xx){
dim <- xx@modes
m <- apply(xx@data, c(2:xx@num_modes), mean)
mm <- aperm(array(m, c(dim[-1],dim[1])), c(xx@num_modes,c(1:(xx@num_modes-1))))
return(xx - mm)
}


## Try the tensorTS package in your browser

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

tensorTS documentation built on Aug. 10, 2021, 9:07 a.m.