CP_TTS: Estimating the tensor time series CP-factor model

Description

CP_TTS() deals with the estimation of the CP-factor model for tensor time series. Let $\mathcal{Y}_t$ be a tensor in $\mathbb{R}^{d_1 \times \cdots \times d_m}$. The tensor CP-factor model is given by $$ \mathcal{Y}_t = \sum_{i=1}^r w_i f_{t,i}\,\mathbf{a}_{i,1} \circ \mathbf{a}_{i,2} \circ \cdots \circ \mathbf{a}_{i,m} + \mathcal{E}_t, \quad t \ge 1, $$ where $1 \le r \le \min_{j \in [m]} d_j$ is a fixed but unknown constant, $\mathcal{E}_t \in \mathbb{R}^{d_1 \times \cdots \times d_m}$ is the idiosyncratic error tensor, $\mathbf{f}_t = (f_{t,1}, \ldots, f_{t,r})'$ is the $r$-dimensional factor vector, and $\mathbf{a}_{i,j}$ is a $d_j$-dimensional loading vector corresponding to the $i$-th factor and the $j$-th mode. Without loss of generality, we assume $ | \mathbf{a}_{i,j} |_2 = 1$ for $i \in [r]$ and $j \in [m]$. This function aims to estimate $r$ and the loading vectors $\{\mathbf{a}_{i,j}\}_{i \in [r], j \in [m]}$ using the method proposed in Chang et al. (2026+).

Usage

CP_TTS(Y, xi = NULL, r = NULL, A.init = NULL, control.DPI = list())

Value

The function returns a list containing the following components:

A.hat: The final iterative loading matrices.
A.init: The initial loading matrices used to start the iteration.
Sigma.yij.xii.1: The thresholded moment vectors used in the iterative updates and inference.
r: The number of factors used in the iterative procedure.
iter.step: The number of iterations performed.
fnorm.resid: The relative Frobenius norm of the residuals recorded during the iterations.
f.hat: The estimated factor series based on the final iterative loading matrices.
f.hat.inl: The estimated factor series based on the initial loading matrices.
delta.sel: The selected threshold level from the initial one-pass estimation. If A.init is supplied by the user, this value is NULL.
control.DPI: The control list actually used in the function after merging user-supplied values with the defaults.

Arguments

Y

An array representing a tensor-valued time series with dimension $n\times d_1 \times \cdots \times d_m$, where $n$ is the sample size and $m \ge 2$.

xi

An auxiliary scalar series $(\xi_1,\ldots,\xi_n)'$, which is a linear combination of $\mathrm{vec}(\mathcal{Y}_t)$. If xi = NULL (the default), $\xi_t$ is estimated by the PCA method described in Chang et al. (2023), or by a randomized projection method by setting random.projection = TRUE in control.DPI.

r

The prescribed number of factors. If set to NULL (the default), $r$ is estimated from the data by the ER method or the Log-ER method by setting ratio.type in control.DPI to "classical" or "log", respectively.

A.init

Optional initial loading matrices. It should be a list of length $m$, where the $j$-th sublist is a $d_j \times r$ matrix. If NULL, an initial estimator is obtained from an initialization step.

control.DPI

A named list of control parameters used in the double projection iteration (DPI) algorithm. The supported components are:

lag.k.dpi: Positive integer. Number of lags $K$ used in $\tilde{\mathbf{M}}_j = \sum_{k = 1}^K \tilde{\mathbf{\Sigma}}_{k,j} \tilde{\mathbf{\Sigma}}_{k,j}'$, where $\tilde{\mathbf{\Sigma}}_{k,j}$ is an estimate of the cross-covariance between $\mathbf{Y}_{t,j}$, the mode-$j$ matricization of $\mathcal{Y}_t$ with dimension $d_j \times \prod_{j' \neq j} d_{j'}$, and $\xi_t$ at lag $k$. Default is $10$.

threshold

Logical. Whether thresholding is applied in the initialization and iteration steps. Default is TRUE.

delta

Optional thresholding level used in the initialization step. Default is NULL. If NULL, it is selected via a grid search method.

delta2

Numeric vector of length m, controlling the thresholding level in each tensor mode during the iterative update. The default j-th element is $\hat{\sigma}_0 (n^{-1}\log d_j)^{1/2}$, where $\hat{\sigma}_0^2 = (n \prod_{j=1}^m d_j)^{-1} \sum_{t=1}^n\|\mathcal{Y}_t\|^2_{\mathrm{F}}$.

ratio.type

Character string specifying the ratio criterion used in estimating $r$. Typical choices are "log" for the Log-ER method and "classical" for the ER method. Default is "log".

random.projection

Logical. If TRUE, a randomized projection step is used to select xi. See Section 3.4 of Chang et al. (2026+) for details. Default is FALSE.

iter.max

Maximum number of iterative updates. Default is $20$.

eps

Stopping tolerance for the iterative algorithm. Default is $10^{-4}$.

grid.num

Integer. Number of grid points used when selecting the thresholding level in the initialization step. Default is $50$.

delta.max

Maximum value of the thresholding grid in the initialization step. Default is $0.1 \hat{\sigma}_0 (n^{-1} \sum_{j = 1}^m \log d_j)^{1/2}$.

print.eps

Logical. Whether to print the iterative convergence measure. Default is FALSE.

iter.lag

Positive integer. Number of candidate lags used in each iterative update. Default is $1$.

all.put

Logical. If TRUE, the iterative routine returns full intermediate outputs; otherwise, only a compact result is returned. Default is FALSE.

A

Optional true loading matrices, used only for diagnostic purposes in simulations. Default is NULL.

component

Optional true common component tensor, used only for diagnostic purposes in simulations. Default is NULL.

Details

The initial method involves the estimation of the autocovariance between $\mathbf{Y}_{t,j}$ and $\xi_t$ at lag $k$, which is defined as follows: $$ \hat{\mathbf{\Sigma}}_{k,j} = T_{\delta_1}\{\hat{\boldsymbol{\Sigma}}_{\mathbf{Y}_j,\xi}(k)\} \quad \mbox{with} \quad \hat{\boldsymbol{\Sigma}}_{\mathbf{Y}_j,\xi}(k) = \frac{1}{n-k}\sum_{t=k+1}^n (\mathbf{Y}_{t,j}-\bar{\mathbf{Y}}_j)(\xi_{t-k}-\bar{\xi}), $$ where $\bar{\mathbf{Y}}_j = n^{-1}\sum_{t=1}^n \mathbf{Y}_{t,j}$, $\bar{\xi} = n^{-1}\sum_{t=1}^n \xi_t$, and $T_{\delta_1}(\cdot)$ is a threshold operator defined as $T_{\delta_1}(\mathbf{W}) = \{w_{i,j}1(|w_{i,j}| \ge \delta_1)\}$ for any matrix $\mathbf{W}=(w_{i,j})$, with threshold level $\delta_1 \ge 0$ and $1(\cdot)$ denoting the indicator function. Chang et al. (2026+) suggest choosing $\delta_1$ by a grid search method. See Section 3.4 of Chang et al. (2026+) for details.

References

Chang, J., He, J., Yang, L., & Yao, Q. (2023). Modelling matrix time series via a tensor CP-decomposition. Journal of the Royal Statistical Society Series B: Statistical Methodology, 85, 127--148. tools:::Rd_expr_doi("10.1093/jrsssb/qkac011").

Chang, J., Huang, G., Yao, Q., & Yu, L. (2026+). CP-factorization for high dimensional tensor time series and double projection iterations. Preprint. tools:::Rd_expr_doi("doi:10.48550/arXiv.2606.08560").

Examples

Run this code

n <- 200
D <- c(10, 10)
r <- 2

data <- HDTSA:::DGP.TCP(
  n = n,
  m = 2,
  D = D,
  r = r,
  w = c(10, 10),
  ar.coef = list(0.5, 0.3),
  factor.loading = "sparse-random",
  alpha = 0.3
)

Y <- data$Y

fit <- CP_TTS(Y)

fit$r
fit$A.hat
fit$f.hat

fit.known <- CP_TTS(Y, r = 2)

Run the code above in your browser using DataLab