R/tenAR.R

Defines functions .remove.mean .getrank .getpos predict.rolling tenAR.predict mplot.acf mplot tenAR.bic tenAR.SE.MLE tenAR.SE.LSE MAR1.CC MAR1.RR tenAR.MLE tenAR.LS tenAR.PROJ tenAR.A tenAR.VAR MAR.SE MAR1.MLE MAR1.LS MAR1.PROJ MAR1.RRCC.SE MAR1.RRLS.SE matAR.RR.se matAR.RR.est tenAR.est tenAR.sim

Documented in matAR.RR.est matAR.RR.se mplot mplot.acf tenAR.est tenAR.predict tenAR.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}.
#' See "Details" of \code{\link{tenAR.est}}.
#'@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.
#'@seealso \code{\link{tenFM.sim}}
#'@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 {
      stop("Please specify cov")
    }
    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 {
    stop("please specify method.")
  }
}


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 n.head prediction horizon.
#'@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){
    return(predict.rolling(A, xx, n.head, method, n0))
  }

  P <- length(A)
  R <- length(A[[1]])
  K <- xx@num_modes - 1
  dim <- xx@modes
  ttt <- (dim[1]+1):(dim[1]+n.head)
  xx.pred <- array(0, c(n.head, prod(dim[-1])))
  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))
  }
  return(array(xx.pred, c(n.head, dim[-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)}
  ttt <- n0:(t - n.head)
  xx.pred <- array(0, c(t-n.head-n0+1, prod(dim[-1])))
  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))
  }
  return(array(xx.pred, c(t-n.head-n0+1, dim[-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.