Threshold Tensor Factor Model in CP Form

Threshold Tensor Factor Model in CP Form

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

} \fi

{\it JEL Classifications}: C13, C32, C55

{\it Keywords:} CP Decompositions, Factor models, Threshold auotoregressive models, High-dimensional, Tensor data.

\spacingset{1.7}

Introduction

Tensor time series data, denoted as $\mathcal{X}_t \in \mathbb{R}^{d_1 \times d_2 \times \cdots \times d_K}$ for $t = 1, \ldots, T$, is a collection of multidimensional arrays observed sequentially over time. Such tensor time series arise in various scientific and financial applications, including multi-category product import-export volume among a group of countries over time ChenYangZhang2022, quarterly economic indices of different countries ChenChenTsay2019, returns of portfolio constructed by size, book to market ratio, and momentum LiXiao2021, multivariate spatial-temporal data chen2020modeling,barigozzi2025general, and many others.

Bolivar2025Review provides a review of recent developments in autoregressive and factor modeling for such data. For example, ChenXiaoYang2021, han2023rr and li2024cointegrated extend the autoregressive model to matrix-valued time series (MAR), while LiXiao2021 and WangZhengLi2024 extend it further to tensor time series. Another widely used approach in high-dimensional data is dimension reduction through factor models. The two main decompositions for tensor data structures are Tucker decomposition and CP decomposition KodaBader2009. Mimicking Tucker decomposition, WhangLiuChen2019, ChenYangZhang2022 and ChenFan2023 introduced a tensor factor model in Tucker form, with factors in tensor form but of much smaller dimensions than the observed tensors. Other notable contributions in this area include ChenChenTsay2019,han2024tensor,han2022rank,chen2025diffusion, and others.

In contrast, han2024cp propose a tensor factor model in CP form (TFM-cp) that uses a set of univeriate factor processes that are uncorrelated, and rank-1 tensor factor loadings. This approach mimics the CP decomposition. This model structure facilitates the analysis of underlying dynamics in high-dimensional tensor time series through a set of univariate uncorrelated latent factor processes, making it significantly more convenient for studying the dynamics of the time-series. The TFM-cp is specified as:

where each $f_{jt}$ is a one-dimensional latent factor, ${{\boldsymbol{u} }}_{jk} \in \mathbb{R}^{d_k}$ are the loading vectors normalized to $\|{{\boldsymbol{u} }}_{jk}\|_2=1$, ${\cal E}_t$ is a zero mean idiosyncratic noise tensor with potentially weak cross-correlations among its elements and assumed to be uncorrelated with the latent factors, and $r$ is the rank of the decomposition, or equivalently the number of factors. Without loss of generality, it is assumed $\mathbb{E} f_{jt}^2\le C<\infty$, for all $1\le j\le r$. Then, the signal strengths are contained in $\lambda_j$. Here, $\{{{\boldsymbol{u} }}_{jk}, 1\le i\le r\}$ are not necessarily orthogonal. chang2023modelling,han2024cp showed that the model can be very useful in many applications.

We observe that the existing approaches mentioned above are all linear models. It has been widely acknowledged that in many applications linear time series models are not sufficient to capture various nonlinear phenominons encountered and there has been a large literature in study nonlinear time series models Tong90, tsay2018nonlinear,Gooijer2017NonLin_book. Modern machine learning algorithms such as deep neural networks are also large nonlinear systems and have shown to be extremely powerful in providing accurate predictions zhou2025factor,luo2025supervised.

Motivated by observed regime-switching phenomenon in applications and encouraged by its simple form, regime-switching models, particularly the threshold models, are a class of powerful nonlinear time series models that have been studied extensively and have been shown to be a simple, elegant and useful parametric model Tong90, Chan&Tong86, Sims&Zha2006, tsay1989testing, ang2012regime,hansen2011threshold in many applications ranging from economics, biology and ecology, signal processing, and many others. As an attempt to move from the linear models to nonlinear ones for tensor time series, in this paper we propose a Threshold dynamic Tensor Factor Model in CP form (T-TFM-cp), which extends the standard TFM-cp model by incorporating a threshold mechanism into the latent factor dynamics. Specifically, we model the factors $f_{jt}$ using a general Threshold Autoregressive (TAR) model:

The TAR model allows the factor dynamics to switch between different regimes depending on the value of a threshold variable $z_{jt}$. This variable can be a function of observable variables (at time $t$) or lagged latent factors. The parameters $p_j$ and $s_j$ determine the specific form of the TAR model for factor $j$.

There is a vast literature on TAR models Tong90, Chan&Tong86, Sims&Zha2006, tsay1989testing, ang2012regime. For matrix autoregressive models, recent work such as Bucci2024, WuChan2024, and YuLiZhangTong2024 develops new methods for modeling regime-specific changes in the coefficient matrices. Related but conceptually distinct, LiuChen2020 introduce a threshold vector factor model, and LiuChen2022 study a threshold matrix factor model in Tucker form. In their framework, the loading matrices differ across regimes and regime switching is controlled by a threshold variable. The latent factors are assumed to be stationary but are not equipped with an explicit dynamic model, and therefore these approaches do not provide forecasting capability. In contrast, the proposed T-TFM-cp model assumes common loading structures across all regimes, while allowing the latent factor processes to switch between regimes according to one or more threshold variables. A distinctive feature of T-TFM-cp is its flexibility: all factors may share the same threshold variable, or each factor can be governed by its own thresholding rule. Moreover, unlike the threshold factor models mentioned above, T-TFM-cp incorporates an explicit dynamic specification for the factor processes, thereby enabling prediction, similar in spirit to the dynamic matrix factor model of YuChen2024.

In this paper, we extend the TFM-cp model to a nonlinear threshold framework--rather than extending the TFM-Tucker model--for several reasons. First, the TFM-cp model (ref) is uniquely defined up to sign changes when the signal strengths are ordered, facilitating interpretation. Second, by decomposing the factor processes into uncorrelated univariate time series, TFM-cp enables a more flexible modeling of temporal dynamics. Compared with Tucker-based representations, CP models typically require fewer factors to capture the essential signal structure, achieving a favorable balance between parsimony and explanatory power. Third, the uncorrelated univariate factor structure enables independent modeling of regime-switching dynamics for each factor, potentially with different threshold variables, threshold values, or numbers of regimes.

By combining CP decomposition with TAR dynamics, the proposed model offers a powerful and flexible framework for analyzing high-dimensional tensor time-series data exhibiting regime-switching behavior in the latent factors. Balancing parsimony, interpretability, and statistical efficiency, it provides a valuable tool for researchers working with complex multiway temporal data in fields such as neuroscience, finance, and signal processing.

The remainder of the paper is organized as follows. Section (ref) presents the proposed model in detail, including the TAR dynamics of the latent factors and the estimation procedures. Section (ref) establishes the theoretical properties of the estimators. Sections (ref) and (ref) demonstrate the effectiveness of the model through simulation studies and an empirical application, respectively. Section (ref) concludes with a discussion of future research directions.

Model Setting and Estimation

In this section, we focus on a more specific specification of the T-TFM-cp model introduced in (ref) and (ref). In particular, we assume that each latent factor follows a threshold autoregressive (TAR) process of the form

where $L$ is the number of regimes, $s_{j,l}$ are the threshold values, $p_j$ are the autoregressive orders, $\phi_{j}^{(l)}$ are the autoregressive coefficients for regime $l$, and $\xi_t$ are white noise processes with variances 1 for $l = 1, 2, \cdots, L$. In the Self-Exciting TAR (SETAR) case, the exogenous threshold variable ${{\boldsymbol{z} }}_t$ is assumed to be past values of $f_{jt}$, e.g. $f_{j,t-\tau_j}$, where $\tau_j$ is the delay parameter. In the following we use the notations ${{\boldsymbol{s} }}_j=(s_{j,1},...,s_{j,L-1})^{^\oldtop} \in \mathbb{R}^{L-1}$, ${{\boldsymbol{\phi} }}_j^{(l)}=(\phi_{j,1}^{(l)},...,\phi_{j,p_j}^{(l)})^{^\oldtop} \in \mathbb{R}^{p_j}$, and ${{\boldsymbol{\phi} }}_j=({{\boldsymbol{\phi} }}_j^{(1){^\oldtop}},...,{{\boldsymbol{\phi} }}_j^{(L){^\oldtop}})^{^\oldtop} \in \mathbb{R}^{L p_j}$.

The TAR model is a powerful and elegant nonlinear extension of the linear autoregressive (AR) model, capable of capturing regime switching driven by an observable threshold variable. It has been widely used to model economic and financial data exhibiting asymmetric dynamics and nonlinear patterns. By allowing different autoregressive behaviors depending on the regime of the threshold variable $z_{j,t-d}$, the model can represent phenomena such as asymmetric cycles, shifts in mean levels, and volatility clustering. This flexibility makes TAR models particularly suitable for characterizing distinct dynamics during economic expansions and recessions Tong90, Chan&Tong86, tsay1989testing.

We refer to the combination of the TFM-cp model in (ref) and the dynamic TAR specification in (ref) as the TAR Tensor Factor Model in CP form (T-TFM-cp). To estimate the parameters of the T-TFM-cp model, we adopt a two-step estimation strategy. In the first step, we estimate the TFM-cp model while ignoring the dynamic structure of the latent factors. This yields estimates of both the loading vectors $\mathbf{u}_{jk}$ and the factor process $\widehat{f}_{jt}$. The procedure follows chen2024estimation, which combines a composite PCA initialization with an iterative simultaneous orthogonalization (ISO) algorithm. Details and theoretical properties can be found in chen2024estimation. Denote the estimated loading vectors by $\hat{{{\boldsymbol{u} }}}_{jk}$ and the estimated factors by $\widehat f_{jt}={\cal X}_t\times_{k=1}^K \hat {{\boldsymbol{v} }}_{jk}^{{^\oldtop}}$, where $\hat{{{\boldsymbol{v} }}}_{jk}$ is the $j$-th column of $\widehat{{{\boldsymbol{U} }}}_k(\widehat{{{\boldsymbol{U} }}}_k^{^\oldtop}\widehat{{{\boldsymbol{U} }}}_k)^{-1}$, with $\widehat{{{\boldsymbol{U} }}}_k=(\hat{{{\boldsymbol{u} }}}_{1k},\cdots,\hat{{{\boldsymbol{u} }}}_{rk})\in \mathbb{R}^{d_k\times r}$. Under this definition, $\widehat f_{jt}$ serves as an estimate of $\lambda_j f_{jt}$ in (ref), implicitly absorbing the signal strength $\lambda_j$. This scaling does not affect the results of the simulation study or the real data analysis; it only changes the magnitude of the estimated factors.

Second, once the factor process estimates $\widehat{f}_{jt}$ are obtained, we fit an individual TAR model to each estimated factor process using standard identification and estimation techniques for TAR models tong1990non,tsay1989testing,tsay2018nonlinear. In particular, given the threshold variable, the number of regimes, and the AR orders of each regime, the AR coefficients and threshold values are estimated by least squares. Specifically, for each $j=1,\ldots,r$, the estimates $\hat{{{\boldsymbol{\phi} }}}_j$ and $\hat{{{\boldsymbol{s} }}}_j,j=1,\ldots,r$ minimizes

Although the estimation of the threshold vector ${{\boldsymbol{s} }}_j$ is strictly a non-convex optimization problem, it is typically not computationally intensive since one only needs to check the observed values of the threshold variable; all values between two consecutive observations of the threshold variable are equivalent chan1993consistency,hansen2000sample. Sequential update of the least squares (when adding or dropping one observation) makes exhaustive search computationally efficient, though more advanced algorithms are available li2016nested,li2011least. The number of regimes and the AR orders are typically determined by BIC or similar measures, as well as diagnostic model checking procedures involving residual analysis, testing for remaining nonlinearity, and evaluating forecasting performance Tong90,hansen2000sample,li2011least,gonzalo2002estimation,chan1998limiting.

The two-step estimation strategy is standard for models consisting of multiple components, such as the T-TFM-cp model. Similar procedures have been widely used in the literature on dynamic (vector) factor models; see, for example, stock2016, Jasiak2001, hallin2016, otto2022approximate.

Threshold variable determination is essential for building a threshold model. When there is no prior knowledge on the threshold variable, a data-driven procedure is needed in order to search for a suitable one. In standard univariate threshold models, a typical candidate pool is the lag variables Tong&Lim1980,Tong90,chan1993consistency,Tong2010 and identification is commonly performed by comparing a handful of plausible specifications. When multiple time series are involved, as in our setting, it may be desirable to identify one or a small number of common threshold variables, which can substantially improve model interpretability. When one considers a large candidate pool consisting of exogenous variables, a trial-and-error approach can be extremely time consuming, complicated more by the multiple comparison problem at the end. In such cases, the reverse approach proposed in Wu&Chen2007,liu2016regime provides an effective alternative. It is also worth noting that different factor processes may rely on different types of threshold variables: some may follow self-exciting dynamics using lagged values of the factor itself, while others may require observable exogenous threshold variables.

Theoretical Properties of the Estimators

In this section, we present some theoretical properties of the proposed two-stage estimation procedures. We introduce some notations first. Let $d=\prod_{k=1}^K d_k ,d_{\max}=\max\{d_1,...,d_K\}$. The matrix Frobenius norm is defined as $\|{{\boldsymbol{A} }}\|_{\rm F} = (\sum_{ij} a_{ij}^2)^{1/2}$. Define the spectral norm as $$ \|{{\boldsymbol{A} }}\|_{2} = \max_{\|{{\boldsymbol{x} }}\|_2=1,\|{{\boldsymbol{y} }}\|_2= 1} \|{{\boldsymbol{x} }}^{{^\oldtop}} {{\boldsymbol{A} }} {{\boldsymbol{y} }}\|_2.$$ Considering that the loading vector ${{\boldsymbol{u} }}_{jk}$ of TFM-cp can only be identified with a change in sign, we employ

to quantify the discrepancy between $\widehat {{\boldsymbol{u} }}_{jk}$ and ${{\boldsymbol{u} }}_{jk}$.

To establish the asymptotic properties of the proposed procedures, we impose the following assumptions.

Assumption (ref) closely resembles the noise conditions found in foundational studies such as bai2002, bai2003, lam2011, LamYao2012, and other notable contributions in the factor model literature. For ease of exposition, we assume that the noise tensor is independent across time $t$, hence allowing for weak cross-sectional contemporaneous dependence. Although introducing a weak temporal correlation among noises, as proposed by bai2002, is feasible, it significantly complicates the theoretical framework. Additionally, we adopt the normality assumption for technical convenience, although it can be generalized to accommodate exponential-type tail conditions, as seen in chen2024estimation.

Assumption (ref) is a standard condition in factor models, such as chen2024estimation and han2024cp. In the case of strong factor models, we have $\lambda= \sqrt{d}$. Assumption (ref) requires the tail probability of $f_{jt}$ to exhibit exponential decay, which is a standard assumption. In particular, when $\gamma_2 = 2$, this implies that $f_{jt}$ follows a sub-Gaussian distribution. The mixing condition is a well-established assumption that encompasses a wide variety of time series models, including causal ARMA processes with continuous innovations, as discussed in Tong90, tsay2005analysis, fan2008nonlinear, rosenblatt2012markov, tsay2018nonlinear, among others.

Next, we impose assumptions for the TAR structure of the factor processes, depending on whether the threshold variable $z_{jt}$ is self-exciting (i.e., the past values of $f_{jt}$) or is an exogenous observable variable. If $z_{jt}$ is self-exciting with $z_{jt}=f_{j,t-\tau_j}$, then ${{\boldsymbol{Y} }}_{jt} = (f_{j,t-1},\ldots,f_{j,t-(p_j\vee\tau_j)})$ is a Markov chain. Denote its $m$-step transition probability by $\mathsf{P}^m (y,A)$, where $y \in \mathbb{R}^{p_j\vee \tau_j}$ and $A$ is a Borel set. If $z_{jt}$ is observable, let ${{\boldsymbol{Y} }}_{jt} = (f_{j,t-1},\ldots,f_{j,t-p_j},z_{jt})$, Denote its $m$-step transition probability by $\mathsf{P}^m (y,A)$, where $y \in \mathbb{R}^{p_j+1}$ and $A$ is a Borel set.

Assumptions (ref), (ref), (ref) are commonly used in theoretical analysis of TAR models. Under Assumption (ref), ${{\boldsymbol{Y} }}_{jt}$ is V-uniformly ergodic with $V(\cdot) = K(1 + \|\cdot\|)$, a condition that is stronger than geometric ergodicity. For a detailed explanation of V-uniform ergodicity, see Chapter 16 in meyn2012markov. In the specific case where all threshold variables are lags of $f_{jt}$ and errors are homoscedastic across all regimes, Assumption (ref) combined with the condition $\max_{\ell} \sum_{i=1}^{p_j} |\phi_{j,i}^{(\ell)}|<1$ suffices for Assumption (ref) to hold. More details can be found in chan1985use and chan1989note. Assumption (ref) implies that the autoregressive function is discontinuous at the thresholds $s_{j,\ell}$. In Assumption (ref)(i), if $\tau_j>p_j$, then $w_{j,\ell,\tau_j}$ may not be a component of ${{\boldsymbol{w} }}_{j\ell}$. In this scenario, Assumption (ref)(i) is equivalent to the conditions ${{\boldsymbol{\phi} }}_j^{(\ell)}\neq {{\boldsymbol{\phi} }}_j^{(\ell+1)}$ for $1\le \ell \le L-1$, which are both necessary and sufficient for the identification of all thresholds. Assumption (ref) (ii) is similar to Assumption 3.6 in zhang2024least.

han2024cp proposed iterative procedure ISO to estimate the factor loading vectors of TFM-cp using auto-covariance tensor. When the number of factors $r$ is fixed, the procedure achieves convergence rate $\| \widehat{{\boldsymbol{u} }}_{jk}\widehat{{\boldsymbol{u} }}_{jk}^{{^\oldtop}} -{{\boldsymbol{u} }}_{jk} {{\boldsymbol{u} }}_{jk}^{{^\oldtop}} \|_{2} =O_{\mathbb{P}}(\sigma d_{\max}^{1/2}\lambda^{-1}T^{-1/2} +\sigma^2\lambda^{-2}),1\le k\le K$. Using their iterative procedure, we can further establish the theoretical properties of the estimated latent factor process $\widehat f_{jt}$ in (ref).

The proposition establishes a non-asymptotic bound for the estimated factors $\widehat f_{jt}$. It demonstrates that consistent factor estimation requires the signal-to-noise ratio ($\lambda/\sigma$) to increase to infinity, ensuring sufficient information about the signal at each time point $t$. As expected, a higher signal-to-noise ratio leads to faster convergence, as stronger factors provide more information in the observed data. Additionally, (ref) shows that the error rate between the sample (auto-)covariance of the estimated factors and the true sample (auto-)cross-moment is $o_{\mathbb{P}}(T^{-1/2})$ when $\lambda/\sigma\gg \sqrt{d_{\max}}+T^{1/4}$. In the case of a strong factor model setting bai2002,lam2011, where $\lambda/\sigma\asymp \sqrt{d}$, this condition is equivalent to $T\ll d^2$, which generally holds. This suggests that under such conditions, using the estimated factor processes for model building and inference of the TAR component is equivalent to using the underlying true factors, without any efficiency loss. The statistical rates in (ref) and (ref) provide a basis for further modeling of the estimated factors.

The following theorems present the convergence rate and asymptotic normality for the least squares estimators (ref) in the second-stage TAR modeling of the latent factor process, using $\widehat f_{jt}$ as the true factor process. Theorems (ref) and (ref) separately examine the cases of self-exciting threshold variables and observable threshold variables.

In the case of self-exciting threshold variables, the estimation error of the threshold level ${{\boldsymbol{s} }}_j$ in (ref) comprises two components: the conventional $1/T$ rate from chan1993consistency,li2012least and zhang2024least and the estimation error of the factor process from the first stage in (ref). This additional error arises because the lagged estimated factor process serves as the threshold variable. Since the estimated factors determine the Markovian regime at each time $t$, we require a non-asymptotic convergence rate for all $\widehat f_{jt}$, as established in Proposition (ref). The convergence rate provided by Theorem 3 in han2024cp is insufficient for this purpose. In contrast, for observable threshold variables, the threshold level ${{\boldsymbol{s} }}_j$ in (ref) achieves the standard $1/T$ rate without additional factor estimation error.

The presence of factor estimation error in the self-exciting threshold variable case leads to larger AR coefficient estimation errors compared to the observable threshold variable case. For observable thresholds, the error consists of the parametric rate $T^{-1/2}$ and the error from the sample covariance matrix using estimated factors as in (ref). For self-exciting thresholds, an additional term $\sigma/\lambda$ appears, reflecting the signal-to-noise ratio. Consequently, the least squares estimator $\widehat {{\boldsymbol{\phi} }}_j $ satisfies a central limit theorem only when the signal strength is sufficiently high: specifically, $\lambda/\sigma\gg \sqrt{d_{\max}}+ T^{1/4}$ for observable threshold variables and $\lambda/\sigma\gg \sqrt{d_{\max}}+ T^{1/2}$ for self-exciting threshold variables. When the signal is weaker, consistency remains achievable, but the factor loading and idiosyncratic noise ${\cal E}_t$ dominate the estimation error. In such regimes, the factor process error in (ref) and/or (ref) dominates the parametric rate, preventing the derivation of a tractable asymptotic distribution.

Simulations

We conducted a simulation study to assess the performance of the proposed T-TFM-cp model in matrix form across a range of settings. The scenarios systematically vary the dimensions, sample sizes, and signal-to-noise ratios (SNR, $\lambda/\sigma$). Specifically, we generate data from weak factor model with $\lambda=1$,

where, for $j=1,2,3$, the components $\boldsymbol{u}_{1,j}$ and $\boldsymbol{u}_{2,j}$ are generated by drawing i.i.d. Gaussian entries with zero mean and then normalizing the resulting vectors. All elements of ${{\boldsymbol{E} }}_t$ have variance $\sigma^2$. The latent factors follow the TAR models

where $a_{j,t}$ are independent Gaussian white noise with zero mean and unit variance ($\varsigma_j^2 = 1$), for $j = 1, 2, 3$. All results are based on 100 simulation replications.

Given that the TFM-cp factors are identifiable only up to sign changes and permutations, we first aligned the estimated factors with the true simulated factors using optimal sign matching and permutation selection. We then evaluate the estimation accuracy of the loading vectors using the Mean Squared Error (MSE):

We use the estimation procedure of han2024cp to estimate the TFM-cp model, using lag 1 ($h=1$) auto-moments in the estimation. Figure (ref) presents boxplots of the logarithm of $\text{MSE}_j$ for the three loading vectors, under three choices of matrix dimensions, three SNR levels, and three time-series lengths. As expected, the estimation error decreases as the sample size, dimensionality, and SNR increase.

Next we study the impact of using the estimated factor process as well as its lag variable as the threshold variable, compared to using the underlying true factor process and its lag threshold variable. To isolate the impact, we fix the true TFM-cp rank $r=3$, use the true AR order each factor as in (ref), and set the threshold delay to its true value ($\tau_j=1$). Treating the estimated $\hat{f}_{j,t}$ ($j=1,2,3$) as observed time series, we estimate the threshold model (the AR parameters and the threshold value) under the true configuration and the estimated threshold variable $\hat{f}_{j,t-1}$. Given the estimated threshold variable $\hat{f}_{j,t-1}$ and the estimated threshold $\hat{s}_{1,j}$, each observation $\hat{f}_{j,t}$ is classified into one of the two regimes (lower vs. upper). We then compute the proportion of observations for which the inferred regime matches the underlying true regime (i.e. $f_{i,t-1}<0$ for the lower regime and $f_{i,t-1}\geq 0$ for the upper regime). For comparison, we repeat the procedure using the “true” factor process $f_{i,t}$ and the “true” threshold variable $f_{i,t-1}$ to estimate the AR coefficients and the threshold $\tilde{s}_{1,j}$. Using $f_{j,t-1}$ and $\tilde{s}_{1,j}$, we compute the corresponding correct-classification proportion. Figure (ref) shows the comparison of these two proportions under the same model settings as Figure (ref). As expected, the results based on the true factor process do not vary across different dimensions or SNR levels, since they are independent of the first-stage TFM-cp estimation. In contrast, the results based on the estimated factor processes approach the performance of the true-factor case as the dimension or SNR increases, reflecting the improved accuracy of the first-stage estimation.

Prediction of ${{\boldsymbol{X} }}_{t+1}$ given ${{\boldsymbol{X} }}_{1},\ldots, {{\boldsymbol{X} }}_{t}$ is carried out by first forecasting the factor process. To assess the prediction performance, we generated 200 additional observations for each simulation replicate beyond the original sample size $T$. Without updating the estimated model parameters, rolling one-step ahead forecasts are produced, and prediction accuracy is assessed using the sum of squared prediction errors \[ \sum_{t=T}^{T+199} ||\widehat{{{\boldsymbol{X} }}}_t(1)-{{\boldsymbol{X} }}_{t+1}||^2_F. \] Here, for $t\geq T$, $\widehat{{{\boldsymbol{X} }}}_t(1)=\sum_{j=1}^3\hat{{{\boldsymbol{u} }}}_{1,j} \hat{{{\boldsymbol{u} }}}_{2,j}'\hat{f}_{j,t}(1)$, where $\hat{f}_{j,t}(1)$ denotes the one-step-ahead forecast from the estimated TAR model. These forecasts are based on the TFM-cp factor estimates $\hat{f}_{1,t}, \hat{f}_{2,t}, \hat{f}_{3,t}$, obtained by projecting ${{\boldsymbol{X} }}_{t}$ onto the space spanned by the estimated loading vectors $\hat{{{\boldsymbol{u} }}}_{1,j}, \hat{{{\boldsymbol{u} }}}_{2,j}$, $j=1,2,3$, which were estimated at $T$. For comparison, prediction errors were also computed using the underlying true factor processes $f_{j,t}$ to build the TAR model and forecast $f_{j,t+1}$, which leads to predictions of ${{\boldsymbol{X} }}_{t+1}$. Figures (ref) presents the box plots of the prediction errors, comparing the prediction using estimated loading vectors and the estimated factor process (marked as “Estimated"), and that using the underlying true loading vectors and the underlying true factor processes (marked as “Real"). The results exhibit a pattern similar to that in Figure (ref): when the dimension of ${{\boldsymbol{X} }}_t$ is large, the prediction performance based on estimated quantities approaches that obtained using the true underlying factors.

To eliminate the impact of the unpredictable noise in the observed ${{\boldsymbol{X} }}_{t+1}$, Figure (ref) shows the performance of \[ \sum_{t=T}^{T+199} ||\hat{{{\boldsymbol{X} }}}_t(1)-{{\boldsymbol{M} }}_{t+1}||^2_F, \] where ${{\boldsymbol{M} }}_t$ is the signal part in (ref). The difference between the two approaches becomes slightly more pronounced, but the overall pattern remains consistent.

These simulation results demonstrate the effectiveness of our proposed approach in accurately estimating the loading vectors, identifying regime switching in latent factors, and predicting future observations in tensor time series data. Across a wide range of model settings, including varying dimensionality, sample size, and signal-to-noise ratios, the method consistently delivers reliable estimation and forecasting performance.

Real data example: Multi-country economic indeces

In this section, we demonstrate the proposed T-TFM-cp model by analyzing the matrix time series containing several economic indices from several countries.

Data

We use the dataset from ChenChenTsay2019, which contains quarterly observations from 1990Q4 to 2016Q4 from the 14 countries including: United States of America (USA), Canada (CAN), New Zealand (NZL), Australia (AUS), Norway (NOR), Ireland (IRL), Denmark (DNK), United Kingdom (GBR), Finland (FIN), Sweden (SWE), France (FRA), Netherlands (NLD), Austria (AUT), Germany (DEU). Our analysis focuses on ten macroeconomic indicators: CPGDFD.d2lnsa, CPGREN.d2lnsa, CPALTT01.d2lnsa, IRLT.dlv, IR3TIB.dlv, PRINTO01.dln, PRMNTO01.dln, LORSGPOR.dln, XTEXVA01.GP, and XTIMVA01.GP.

Table (ref) summarizes each series, including its short name, its mnemonic (the series label used in the OECD database), the transformation applied to the series, and a brief description of the data. All data are obtained from the OECD Database. In the transformation column, $\Delta$ denotes the first difference, and $\Delta \ln$ denotes the first difference of the logarithm. GP denotes the growth rate measure from the last period.

Figure (ref) shows the standardized time series data; the rows and the indices by the columns represent countries.

Exploratory Analysis

As an exploratory analysis, univariate ARMA model is fitted to each individual time series, using corrected Akaike Information Criterion (AICc) approach hurvich1989regression,hyndman2008automatic for model determination.

Among the 130 selected univariate ARMA models, 28 are MA(1), 23 are AR(1), 16 are AR of order 4 or 5, 15 are white noise, 10 are ARMA(1,1), and some other less frequently used models. Figure (ref) displays the orders of the ARMA models for each time series selected by AICc. Notably, the first three indices (CPI related) show a different pattern from the remaining seven indices.

We stack each observed matrix into a 130 dimensional vector, and estimate a vector factor model (VFM) with two factors using the estimation procedure of LamYao2012, with $h=1$. The number of factors is determined using the eigen-ratio criterion in LamYao2012. A TFM-cp model is estimated with two factors, again using the estimation procedure of han2024cp with lag 1 ($h=1$) auto-moments. Two factors are used so to be comparable with the VFM model. Note that the loading matrix in the VFM model uses $130\times 2-2$ parameters, while that in TFM-cp uses $(13-1+10-1)\times 2$ parameters for the two (standardized) rank-one $13\times 10$ matrices.

Figure (ref) shows the estimated factors of VFM (left) and TFM-cp (right). The two figures are relatively similar, except for the financial crisis period around the beginning of 2009. During the period, VFM tries to capture the extreme event with one factor, while TFM-cp uses both factors with less magnitude.

To further compare the estimated VFM and TFM-cp, each column of the estimated loading matrix of the VFM model is rearranged into a matrix according to the column and row classifications of the corresponding observed element in the stacked vector. Figure (ref) shows the rearranged matrix of the VFM loading vectors (left) and the estimated rank-one loading matrices of TFM-cp model (right). There are some similarity in both models. For example, the first three indices (CPI related) are loaded heavily on the second factor in both models. Index 4 (long interest rate) is very weakly related to Factor 1, but (negatively) loads on the second Factor in VFM, but not for TFM-cp. Most of the indeces have negative loadings on Factor 2 for VFM, while the CPI related indices have negative loading on Factor 1 for TFM-cp. The TFM-cp model show more country-wide differences in the loading. For example, in TFM-cp, DNK, NOR and AUS are very weakly loaded on Factor 1 in TFM-cp, and NOR is weakly loaded on Factor 2.

For each model, we compute the in-sample fitted values and obtain residuals for each time series. Figure (ref) reports the in-sample coefficient of determination, $R^2=1-\sum(y_t-\hat{y}_t)^2/\sum y^2_t$ for each series under the best ARMA model selected by AICc shown in Figure (ref), as well as the VFM and TFM-cp specifications. The overall in-sample $R^2$ values for the ARMA, VFM, and TFM-cp models are 0.2993, 0.3839, and 0.2937, respectively, using 258, $257 = 130\times 2 - 3$, and $42 = 2{(10-1)+(13-1)}$ parameters.

Figure (ref) indicates that the ARMA models provide markedly better fits for the first three CPI related series—albeit with relatively high ARMA orders (see Figure (ref))—whereas the factor models yield better fits for indices 5–10. None of the models captures the dynamics of the fourth index (IR: Long) satisfactorily.

More specifically, for the first three series, the overall $R^2$ values are 0.5941 (ARMA), 0.3331 (VFM), and 0.3437 (TFM-cp); for the remaining indices, the corresponding values are 0.1633, 0.3470, and 0.3521. It is worth noting that VFM employs substantially more parameters than TFM-cp.

Analysis using T-TFM-cp

The VFM and TFM-cp do not specify any temporal structure of the factors. In the following we present results of fitting ARMA and TAR models for the estimated factors under each setting. The results show that there is indeed a threshold phenomenon among the factors.

Table (ref) reports the models estimated using the first 65 observations; the remaining 40 observations are reserved for evaluating out-of-sample predictive performance. Based on preliminary exploration, we use the log growth of U.S. GDP, $z_t= ln(US.GDP_t) - ln(US.GDP_{t-1})$ as the threshold variable for all factors. U.S. GDP growth is often used as the threshold variable in economic studies as it is a good representative for the status of the economy enders2007threshold,osinska2020modeling,hansen2011threshold,tiao1994some. Over the estimation period, the resulting threshold variable has a mean of 0.0130 and a standard deviation of 0.0050

The T-TFM-cp estimation results reveal markedly different dynamics across regimes for both factors. Factor 1 exhibits a clear threshold effect at 0.0113 (with delay of $d=4$ quarters), splitting the sample into a relatively small regime with 17 observations and a large regime with 44 observations. The upper regime (Regime 1) shows pronounced dynamics, with a strong $AR(1)$ effect ($0.90$) and a significant negative AR(2) term ($-0.43$), indicating substantial oscillatory behavior. In contrast, Regime 2 is well-described by a simple AR(1) structure with a more moderate coefficient ($0.65$), suggesting a smoother dynamics. The threshold variable used for factor 2 is the same as factor 1, with a smaller threshold value, resulting that the upper regime has more observations. Note that the overall residual variances and AIC values of both TAR model for factor 1 and factor 2 are much smaller that of fitting the linear ARMA models to the (same) factors under TFM-cp. The overall variances and AIC values under T-TFM-cp is also smaller than that under VFM using both ARMA and TAR model for the factor processes.

The estimated models using the full sample size is shown in Table (ref) in Appendix Appendix B. The T-TFM-cp model under full sample size is very similar to that of using the first 65 observations. Even though the threshold values are slightly different, the proportions of observations in each regime are similar, showing the relative stability of the estimated models.

Figure (ref) illustrates the regimes identified by the T–TFM-cp model. Red dots indicate Regime 1, and blue triangles indicate Regime 2. The top panel shows the threshold variable for reference, while the middle and bottom panels display Factors 1 and 2, respectively. The two factors exhibit slightly different thresholds, with values of 0.0113 for Factor 1 and 0.0121 for Factor 2, highlighting subtle differences in their regime separation.

To assess the performance of each model listed in Table (ref), out-sample rolling 1-step prediction errors are obtained for the last 40 observations in the data set, from first quarter of 2007 to forth quarter of 2016. In the rolling prediction, the structure of the models (e.g. number of factors, the estimated loading vectors, the ARMA orders or the Threshold AR orders, and threshold variable, the number of regimes and the threshold values) are fixed as listed in Table (ref). Other model parameters, including the factors in the prediction period, the parameters in ARMA or the TAR, are updated using observations 1 to $t$ for the prediction of $f_{j,t+1}$ and ${{\boldsymbol{X} }}_{t+1}$.

{Table (ref)} shows the prediction MSEs over the prediction period under five different models (individual ARMA, ARMA-VFM, ARMA-TFM-cp, T-VFM and T-TFM-cp). In addition, the table also reports the overall MS of the original data in the prediction period, and the in-sample MSE of TFM-cp in the same period, computed using the loading vectors estimated from the full dataset, as references. In the last two rows, we report the prediction performance of the first three indices and the others separately.

From the table, it is seen that T-TFM-cp performed the best considering all indices and indices 4-10 only. The individual ARMA models perform the best for the three CPI indices, similar to that has been revealed in Figure (ref). It is also evident that there is indeed a threshold phenomenon and using the threshold model is useful, as seen from the comparison between ARMA-TFM-cp and T-TFM-cp, and the comparison between ARMA-VFM and T-VFM.

Figure (ref) shows the prediction MSE of each of the 40 prediction times. It seems that for the economic indices 4 to 10, the five models are comparable except during the financial crisis period (mid 2008 to mid 2009), in which T-TFM-cp out-performs ARMA and ARMA-TFM-cp in 2008, and ARMA-VFM in 2009. For the CPI related indices, individual ARMA models perform the best except for the second half of year 2008,

Conclusions

In this paper we proposed a threshold Tensor factor model in CP form, in which the latent factor processes in a TFM-cp model are assumed to follow threshold AR models. It effectively provides a threshold type dynamics for high dimensional tensor time series that allows regime change according to the threshold variable. The model also allows prediction capability while enjoys significant dimension reduction using the factor model structure.

While we have shown that the proposed model is useful in empirical studies, there are challenging issues that are worth further investigation. For example, the two-step estimation procedure works only when the signal to noise ratio in the factor model part is sufficiently high. An effective joint estimation procedure may be needed when the signal to noise ratio is relatively smaller. Hunting for an effective threshold variable is always a challenge in threshold modeling.

\phantomsection