{'fraction_non_alphanumeric': 0.040780770139381056, 'fraction_numerical': 0.013819985825655563, 'mean_word_length': 4.391657379398185, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic Tensor Recommender Systems

Recommender systems have been extensively used by the entertainment industry, business marketing and the biomedical industry. In addition to its capacity of providing preference-based recommendations as an unsupervised learning methodology, it has been also proven useful in sales forecasting, product introduction and other production related businesses. Since some consumers and companies need a recommendation or prediction for future budget, labor and supply chain coordination, dynamic recommender systems for precise forecasting have become extremely necessary. In this article, we propose a new recommendation method, namely the dynamic tensor recommender system (DTRS), which aims particularly at forecasting future recommendation. The proposed method utilizes a tensor-valued function of time to integrate time and contextual information, and creates a time-varying coefficient model for temporal tensor factorization through a polynomial spline approximation. Major advantages of the proposed method include competitive future recommendation predictions and effective prediction interval estimations. In theory, we establish the convergence rate of the proposed tensor factorization and asymptotic normality of the spline coefficient estimator. The proposed method is applied to simulations and IRI marketing data. Numerical studies demonstrate that the proposed method outperforms existing methods in terms of future time forecasting.


Introduction
Recommender systems (RS) are widely used in our daily lives, such as for selecting movies, restaurants, news articles, or online shopping. As one of the information filtering techniques, RS can help users to find interesting items through combining several information sources, e.g., users' ratings and purchasing histories, item profiles and sales volumes, time, location, and companion or promotion strategies. Particularly, incorporating time is useful in RS since users' purchase behaviors are dynamic and often highly dependent on seasonal and time factors, and business sectors also rely on dynamic recommendations to track users' changing purchase interests over time. Thus, it is essential to capture information related to time and develop time-dependent RS, and we refer this as dynamic RS (DRS).
However, developing competitive DRS brings new challenges. First, since data are streaming in over time and are time-dependent, general RS methods which are not capable of capturing time-dependency features may have reduced recommendation accuracy. Second, forecasting future recommendations accurately is also a great challenge for DRS due to the complexity of changing users' interests. For example, users might like to watch news on weekdays, but watch movies on weekends. A shoe store sells more sandals in summer and more snow boots in winter. It is important to borrow information from historical data in developing trends. Many RS methods are not designed to capture trends and predict future recommendations. In addition, as data are streaming in over time, future recommendations could involve new users or new items, whose information is not available from historical data. This is also a common problem encountered in RS, referred as the "cold start" problem.
General RS approaches include content-based filtering and collaborative filtering (CF).
Traditionally, content-based filtering methods recommend similar types of items by matching a user's preferred item profile with current item's profile (e.g., Salter and Antonopoulos, 2006;Son and Kim, 2017). In contrast, CF methods recommend items by predicting item ratings for the active user based on ratings from other similar users (e.g., Herlocker et al., 2004;Luo et al., 2012). On the basis of CF methods, research work related to DRS have been developed in recent years (e.g., Koren, 2009;Gultekin and Paisley, 2014;Yu et al., 2016;Wu et al., 2017;Guo et al., 2018;Xiong et al., 2010;Rafailidis and Nanopou-los, 2014;Bi et al., 2018;Wu et al., 2019). However, most of these methods can only make recommendations for observed discrete time points, and are not designed for future recommendation prediction on unobserved time points; for example, matrix factorization incorporating periodic and continual temporal effects (Guo et al., 2018), coupled tensor factorization exploiting users' demographic information (Rafailidis and Nanopoulos, 2014) and the collaborative Kalman filter (Gultekin and Paisley, 2014). CF methods incorporating a time series model (Yu et al., 2016) or incorporating long short-term memory modeling (Wu et al., 2017(Wu et al., , 2019 are able to solve the forecasting problem, but cannot deal with new users, items or contextual variables. Xiong et al. (2010) used a Bayesian estimation procedure with a time-dependent constraint to predict DRS for new users and items, while Bi et al. (2018) created an additional layer of nested latent factors for new time points, users and items. However, both methods require discrete time points for constructing a tensor.
In this article, we propose a new time-varying coefficient model for the DRS based on tensor canonical polyadic decomposition (CPD); namely, the dynamic tensor recommender system (DTRS). Specifically, we introduce a tensor-valued function of time with each mode corresponding to user, item or a contextual variable, where each component of the tensor is a function of time and has intra-cluster correlation. In the CPD framework, we build a time-varying coefficient model incorporating group information of time points, users, items and contexts. We approximate each coefficient function by a polynomial spline and employ group factors to explore homogeneous group effects. We adopt the weighted least square approach to incorporate intra-cluster correlation for more efficient estimation. In addition, we construct the prediction intervals of estimators of tensor components to forecast the confidence range of predicted values. In theory, we establish the convergence rate of the proposed tensor factorization and the asymptotic property of the spline parametric estimator.
The proposed method has two significant contributions. First, it can effectively forecast recommendations at future time points. This is because the proposed model integrates time dependency to the DRS through time-varying coefficient modeling in tensor factorization so that it can effectively capture dynamic trends of DRS. In addition, the subgroup factors in the proposed model extract homogeneous information from the same group, which provides recommendation forecasting for future time points and therefore solves the "cold start" problem. In contrast to general CF methods which require discrete time points as a tensor mode, the proposed approach is more flexible by utilizing a continuous tensor-value function.
Second, the proposed method is able to provide pointwise prediction intervals. In practice, it is desirable to know the upper and lower bound for predictions, such as the highest possible cost, or the future sales volumes or revenues in the worst case scenario. However, existing methods on prediction intervals are mostly univariate or multivariate time series, and the prediction intervals for user-item-context interactions in a tensor framework have not been developed. The proposed approach develops the prediction intervals for component estimators of a tensor-valued function, which provide a more complete picture of the DRS over time. In our real data analysis, the proposed approach provides effective prediction interval estimators of the sales volumes for IRI marketing data (Bronnenberg et al., 2008), which can help store managers to make sound decisions on marketing strategy and inventory planning.
The remainder of the paper is organized as follows. Section 2 introduces the notation and background on tensor and tensor factorization. Section 3 presents the proposed method and its implementation. Theoretical properties are derived in Section 4. Section 5 presents simulation studies to assess the performance of the proposed approach. In Section 6, we apply the proposed method to the IRI marketing data. Concluding remarks and discussion are provided in Section 7.

Notation and Background
In this section, we introduce the background of the tensor and some notation. Throughout this article, we use blackboard capital letters for sets, e.g., T, I, small letters for scalars, e.g., x, y ∈ R, bold small letters for vectors, e.g., x, y ∈ R n , bold capital letters for matrices, e.g., X, Y ∈ R n 1 ×n 2 , and Euler script fonts for tensors, e.g., X , Y ∈ R n 1 ×n 2 ×···×n d (d > 2).
A dth-order tensor is an array with d dimensions (d > 2), which is an extension of a matrix to higher order. Here d represents the tensor's order. We denote the component (i 1 , i 2 , · · · , i d ) of a dth-order tensor Y by y i 1 i 2 ···i d , where i k = 1, 2, . . . , n k , and k is called a mode of the tensor (k = 1, 2, . . . , d ). In particular, a tensor Y is called a rank-one tensor if it can be written as Y = p 1 • p 2 • · · · p d , where the symbol • represents the vector outer product, and p k is a n k -dimensional latent factor corresponding to the kth mode. That is, each component of the tensor is the product of the corresponding vector components: The canonical polyadic decomposition (CPD) is commonly adopted in tensor decomposition, which decomposes a tensor as a sum of r rank-one tensors. That is: where p k ·j = (p k 1j , · · · , p k n k j ) is a n k -dimensional latent factor corresponding to the kth mode for k = 1, . . . , d; j = 1, . . . , r. Equivalently, each component of Y is The CPD can be considered to be a higher-order generalization of matrix factorisation. Figure 1 illustrates a matrix factorization of a matrix and a CPD of a third-order tensor.
An extensive review of tensors and other forms of tensor decomposition are discussed in Kolda and Bader (2009). into r rank-1 matrices, (b) CPD of a third-order tensor into r rank-1 tensors.
Let P k = (p k ·1 , p k ·2 , . . . , p k ·r ) n k ×r and θ = {P 1 , P 2 , · · · , P d }. We can estimate θ via minimizing a loss function (e.g., L 2 loss). However, the non-convexity of the loss function could impose computational complexity due to numerical instability or even non-convergence (de Silva and Lim, 2008;Frolov and Oseledets, 2017). A common approach to alleviate the non-convexity problem is to introduce regularization. We define an objective function with a penalty function as the following: where Q is a loss function and J is a penalty function, such as L 2 , L 1 or L 0 penalties, or a fused Lasso.
Specially, the optimization problem solves θ * = arg min L(θ|Y), where θ * defines an optimal set of model parameters. In the case of squared loss function with an L 2 -penalty, the objective function is where · F represents the Frobenius norm, and Ω = {(i 1 , i 2 , . . . , i d ) : is a set of indices corresponding to the observed components. Notice that, in the context of RS, the set Ω may not contain all indices of the tensor components and could be a small fraction of the entire tensor size, since the majority of the tensor components could be missing. Major algorithms for implementing the optimization problem include the cyclic coordinate descent algorithm, the stochastic gradient descent method and the maximum block improvement algorithm (Chen et al., 2012).

General Methodology
In this subsection, we develop the methodology for the proposed DTRS method. Specifi-  We consider a dth-order tensor-valued function Y(t) ∈ R n 1 ×n 1 ×...×n d , where the value is the corresponding stochastic process defined on a compact interval T. Without loss of generality, let T be a closed interval [0,1]. Notice that we do not require t i 's to be equally distanced as in other time-dependent tensor models (e.g., Xiong et al., 2010;Bi et al., 2018), so that the proposed method can be applied for a recommender system with arbitrary time points. Figure 2 illustrates an example of a tensor-valued process with d = 3. In the DRS, the tensor-valued process could be the rating or sale volume of items or products from users or stores given contexts. We assume that time points can be categorized into different subgroups, where time points of the same group have common information. For example, in our numerical studies, time points in the same month from the twelve months of each year are categorized in the same group. In addition to time, we also categorize subjects from other modes into subgroups if they share similar characteristics, for example, stores of the same market and products of the same product category.
Suppose the subgroup labels are given, we formulate each component of Y(t) as follows: where ε i 1 i 2 ...i d (t) is a stochastic process with mean zero and finite variance, h j (t) is a trend function of time, p k i k j and q k i k are the jth latent factor and the subgroup factor for the i k th subject from the kth mode, respectively, k = 1, 2, . . . d, j = 1, 2, . . . , r, and g(t) = m d+1 e=1 g e (t)I(t ∈ s e ), in which I(·) is an indicator function, m d+1 is the number of subgroups for time, and g e (t) is a trend function corresponding to the eth subgroup s e of time. We have if the i k th and i k th subjects are from the e k th subgroup (e k = 1, 2, . . . , m k ), where q k (e k ) is the subgroup factor associated with the e k th subgroup, and m k is the number of subgroups for the kth mode, k = 1, 2, . . . , d. We denote the set of observed time points for the component y i 1 i 2 ···i d (t) by T i 1 i 2 ...i d , and the number of components of this set by |T i 1 i 2 ...i d |.
vectors. We assume that the covariance matrix is cov(y typically not an identity matrix due to the intra-cluster correlation arising from repeated observed data. Model (1) adopts the idea of varying-coefficient models to create a CPD for tensor data.
Varying-coefficient models are a useful tool to explore dynamic patterns, and have been applied to modeling and predicting longitudinal, functional, and time series data (Huang and Shen, 2004;Fan and Zhang, 2008). The first part of equation (1) is an individual-level factor model which takes into account the heterogeneity of subjects and trend of time, and the time-varying coefficients h j (t) (j = 1, . . . , r) reflect the dynamic features. The second part of equation (1) is a subgroup-level factor model to capture common features from the same subgroups, where the subgroup factors can accommodate new subjects from any mode at future time points, and the g(t) allows time variables to follow a subgroup function of time such that we can predict future time points via borrowing information from existing time points of the same group.
In addition, since the spline approximation is computationally fast (Xue and Yang, 2006), the model (2) can achieve the spline estimates of the coefficients efficiently, and this is especially advantageous in estimating high-dimensional parameters in RS.
Due to the intra-cluster correlation, it is important to incorporate intra-cluster correlation into RS. However, in practice, the covariance matrix Σ 0 i 1 i 2 ...i d is often unknown. We adopt an invertible working covariance matrix, denoted as (1) , . . . , q k (m k ) ) , and k = 1, . . . , d. Define θ = {P, q, γ} as parameters of interest. We define the following weighted penalized objective function: where λ is the penalized parameter, Ω = {(i 1 , i 2 , . . . , i d ) : y i 1 i 2 ···i d (t) is observed at some t}, N = |Ω| is the number of components of the set Ω, · 2 is the Euclidean norm, and commonly used working correlation structures include independence, exchangeable, and first-order autoregressive process (AR-1), among others. Given a working correlation structure, the working correlation matrix depends on fewer nuisance parameters which can be estimated by the residual-based moment method (Liang and Zeger, 1986). The proposed method is robust to the misspecification of correlation structure as indicated by our numerical examples.

Parameter Estimation
In this subsection, we discuss parameter estimation by minimizing (3) is observed at some t given i k } be the set of indices with the fixed kth mode index i k , where the corresponding components are observed at some time points. We assume that the number of observations for each time subgroup s e is larger or equal than 2 for e = 1, . . . , m d+1 , and the number of observations for each subgroup e k from the kth mode is larger or equal than 2 for e k = 1, . . . , m k ; k = 1, . . . , d. The partial derivatives of L(·|Y) have explicit forms with respect to the individual factors, the subgroup factors and the spline coefficients, which makes it feasible to apply the blockwise coordinate descent approach (BCD). That is, for i k = 1, . . . , n k and k = 1, . . . , d, In fact, the estimation procedure of p k i k in (4) is a ridge regression, and does not require knowing p k i k for i k = i k . Thus, parallel computation is applicable to calculate p k 1 , . . . , p k n k efficiently. The minimization of L(θ|Y) can be done cyclically through estimating P, q, α and β. Notice that Ω = ∪ n k i k =1 Ω k i k , and it is possible that Ω k i k is empty for certain i k 's, that is, there is no observation on the subject i k . Under this circumstance, the individual factor of the i k subject is assigned as p k i k = 0 0, and the predicted values may degenerate to the subgroup-level factor model by utilizing information from members of the same subgroup.

Implementation
In the following, we discuss several implementation issues. To solve the objective function (3), we incorporate the maximum block improvement (MBI) strategy (Chen et al., 2012) into the BCD algorithm cyclically as in Bi et al. (2018). The MBI has two advantages over traditional cyclic BCD algorithms. First, it has a good algorithmic property which guarantees convergence to a stationary point, whereas traditional BCDs may end up with certain points where the criterion function ceases to decrease (Chen et al., 2012). Second, the MBI has the capability of choosing descending directions and hence has the possibility to discover "shortcuts", which may reduce the computational time significantly. Let θ l be an estimator of θ at the lth iteration, θ a be a subset of θ, θ c be the complementary set of θ a , and θ * a be the attempted update of θ a . The improvement of the θ * a is defined as We summarize the implementation of the specifical algorithm as follows.
(i) For each P k , solve (4) through parallel computing and obtain P k * . Then calculate J P k * through (8).
Set the final estimator θ = θ l . Otherwise set l ← l + 1 and go to step 2.
To select tuning parameter λ, we search the one from grid points minimizing the root mean square error on the validation set, defined as where Γ is the set of indices and times of observed data. We choose the number of individual latent factors r such that it is sufficiently large and leads to stable estimation. In general, the r is no smaller than the theoretical rank of the tensor in order to represent subjects' latent features sufficiently well, but not so large as to over-burden the computational cost.
An appropriate selection of the knot sequence is important to efficiently implement the proposed method. In practice, knot locations are usually chosen to be equally-spaced over the range of data or placed at evenly-spaced quantiles of data. Since there are highdimensional factor parameters, for simplicity we set the number of knots to be the integer part of N 1/(2κ+3) , where N = |Ω| and κ is the degree of polynomials. One can also choose other methods to select the number of knots such as the AIC or BIC procedures (Xue and Yang, 2006). The degree of polynomials κ is commonly chosen as 1, 2, or 3. In our numerical study, we set κ = 2 and adopt truncated polynomial bases. One can also use different degrees and spline bases for different time-varying coefficients.
Another important issue is the selection of contextual variables as tensor modes. On the one hand, a higher-order tensor with more contextual variables allows higher-order interactions and hence provides more accurate estimation. On the other hand, a higherorder tensor entails more complex and intensive computation, and may lead to overfitting.
Thus, it is still an important open problem to determine which contextual variables should be included in the tensor. In our numerical studies, promotion strategies are incorporated as a contextual variable, since users' and items' behaviors are distinctive under different promotion strategies. In general practice, however, we assume that the order of a tensor can be determined based on prior knowledge.

Theoretical Properties
In this section, we derive asymptotic properties for the proposed method. Specifically, we establish the convergence rate of the proposed tensor factorization and the asymptotic normality of the spline coefficient estimator. Note that identifiability is critical for tensor representation. We first present the sufficient conditions to ensure identifiability of the proposed tensor modeling as follows.
Proposition 1 If d k=1 K k 2r + d + 1 holds, minimizers of L(P, q, α, β|Y) in P, q, α and β given fixed spline bases are unique up to permutation almost surely, where K k is the Kruskal rank of (P k , q k ), and q k = (q k 1 , q k 2 , · · · , q k n k ) .
Proposition 1 shows that the proposed tensor modeling is identifiable up to permutation almost surely. To address permutation indeterminacy, we could align the factors according to a descending order of the first row of mode-1 factor matrix P 1 , that is, p 1 11 p 1 12 · · · p 1 1r , following the method in Zhang et al. (2014). The rearrangement can be implemented during or after the proposed algorithm, since it does not affect the estimation procedure.
In the rest of Section 4, we assume that the parameters are identifiable.
Thus, the corresponding vector form is Let J(U) be a non-negative penalty function of U. The overall criterion given h j (·) and g(·) is redefined as for U ∈ S, where S is the parameter space for U.
(C1) The functions h j (·) and g e (·) are ξth-order continuously differential for some ξ 2, all j = 1, . . . , d, and e = 1, . . . , m d+1 . The density function of design points t is absolutely continuous and bounded away from zero and infinity on a compact support T.
Conditions (C1)-(C3) are standard in the polynomial spline framework. Similar conditions are also presented in Huang (2003) and Claeskens et al. (2009). In particular, condition (C1) imposes a smoothness condition of trend functions and a mild condition on time density, and guarantees that the observation time points are randomly scattered. Condition (C2) indicates that the adjacent distances among the knot sequence are comparable.

Condition (C3) implies that the eigenvalues of random errors are bounded. Condition (C4)
implies that the difference between the working covariance and true covariance matrices is bounded. Condition (C5) implies that the number of the observed time points grows as the number of the observed components of the tensor increases, to ensure the convergence of the proposed tensor factorization. The following theorem establishes the convergence rate for the proposed tensor factorization.
Theorem 1 Under conditions (C1)-(C5), if the penalty function J(U) has bounded first and second derivatives at true parameter U 0 , as N → ∞, on a δ-ball centered at U 0 for some δ > 0, there exists a minimizer U of (9) such that Theorem 1 provides the convergence rate of the proposed method given trend functions.
When τ = υ, that is, T max and T min have the same order, the convergence rate of the estimator U reaches the optimal rate N −1/2 . Meanwhile, if the order of T max is √ N faster than that of T min , that is, τ − υ = 0.5, then U will not converge to the true U 0 . This implies that to guarantee consistency of the tensor factorization, one should collect sufficient observations even for the least popular user-item-context combinations. In the following theorem, we establish the asymptotic property of the spline coefficient estimator.
Theorem 2 Under conditions (C1)-(C5), if lim N →∞ a N log a N /N = 0 and lim N →∞ a −ξ N N τ = 0, then for any vector c whose components are not all zero, the parametric estimator γ by (6) and (7) satisfies Theorem 2 establishes the asymptotic normality of the spline coefficient estimator. N 1+υ ), and similar results can be found in . The asymptotic variance in Theorem 2 depends on the working covariance matrix and the true covariance matrix. When the working covariance matrices are equal to the true covariance matrices, the asymptotic variance of the proposed estimator reaches the minimum in the sense of Loewner order and the proposed estimator is asymptotic efficient.

The convergence rate of the spline coefficient estimator is
More importantly, the result of Theorem 2 is the key foundation for constructing prediction intervals. First, we derive the standard error for the spline parametric estimates given a fixed λ using the sandwich covariance formula Cov( γ) = ( Ψ + λI) −1 Φ( Ψ + λI) −1 , where , ⊗ operation is the vector operation a ⊗2 = aa , and I is an identity ma- (Chatfield, 1993) where φ σ/2 is the 100(1 − σ)th percentile of the standard normal distribution, and the var{e i 1 i 2 ...i d (t)} is the variance of the prediction error and can be estimated as: The first term in equation (11) is due to estimation error, and the second term can be estimated by the mean squared error on training data.

Simulation Studies
In this section, we perform simulation studies to compare the proposed method (DTRS) with two competing methods, including Bayesian probabilistic tensor factorization (BPTF, Xiong et al., 2010) and the recommendation engine of multilayers (REM, Bi et al., 2018).
The error ε i 1 i 2 i 3 = (ε i 1 i 2 i 3 (t 1 ), . . . , ε i 1 i 2 i 3 (t T )) follows a multivariate normal distribution with mean 0 and a common marginal variance 1, and the correlation structure is either independence or AR-1 with correlation ρ = 0.85.
In each simulation, we consider the number of time points as T = T 1 + T 2 , where the tensor data in the first T 1 = 12 time points are set as the training data, and the tensor data in the last T 2 time points are used as the testing data. For evaluating the forecasting performance at future time points, we consider T 2 = 8 or 12. Considering the missing case, we generate n 1 n 2 n 3 T (1 − π m ) components out of the tensor functions, where π m is the missing percentage and set as 80%. Furthermore, we use π cs = 30% to represent the proportion of new items in the testing data unavailable from the training set. To illustrate the effect of incorporating intra-cluster correlation on estimation efficiency, we compare the estimation efficiency of the proposed methods using different working correlation structures: independent or AR-1, denoted as DTRSin and DTRSar, respectively.
According to Xiong et al. (2010) and Bi et al. (2018), BPTF and REM methods model fourth-order tensor with user, context, item and time modes. For all methods, we assume that the subgroup structure and the number of latent factors are known. For REM and the proposed methods, the tuning parameter λ is pre-selected from grid points ranging from 0 to 20. The validation set is the data from the last four time points of the training set.
For BPTF, we keep the remaining parameters by their default choices. All methods are replicated by 100 simulation runs. Table 1 provides the estimation results of all methods. We observe that the proposed method has better performance when the working correlation structure is the same as the true correlation structure. When the true correlation structure is independence, the DTRSin has smaller RMSE and MAE than the DTRSar, with more than 2.17% improvement. Similarly, when the true correlation structure is AR-1, the DTRSar outperforms the DTRSin. Moreover, the PICPs of the DTRS method are close to 0.95, which implies that the proposed method provides accurate prediction intervals for estimators. For the performance of forecasting time points further away, we observe that the DTRSin, DTRSar and REM methods are relatively robust against time. However, the RMSE and MAE of the REM are larger than the RMSE and MAE of the DTRS method, and the BPTF performs worse on forecasting time points further away. Specifically, the DTRS method performs the best across all settings. For the DTRS method, the relative increasing ratios of RMSEs when T 2 = 8 to those when T 2 = 12 are less than 6.9%, and the corresponding MAEs are at most 4.0%. However, for the BPTF method, the relative increasing ratios of RMSEs for the two time points are more than 9.5%, and the corresponding MAEs are at least 8.1%.
The DTRS method improves on the RMSE and MAE of the BPTF by more than 65%, and improves the RMSE and MAE of the REM by more than 44%. This indicates that the proposed method can obtain more accurate forecasting compared to the BPTF and REM.

Empirical Examples for IRI Marketing Data
In this section, we focus on sales data from drug stores from the IRI Marketing Data (Bronnenberg et al., 2008) to illustrate the performance of the proposed method. The original IRI data is an immense collection of consumer panel data and store sales at grocery stores, drug stores and mass-market stores over the years 2001-2011. The store sales data contain weekly product sales volumes, pricing, and promotion data for all items from 31 product categories sold in 50 U.S. markets. These markets are geographic units defined typically as an agglomeration of counties, usually covering a major metropolitan areas (e.g., Chicago, IL) but sometimes covering just part of a region (e.g., New England). A detailed description of an early version of the data is available in Bronnenberg et al. (2008). types of merchandise display, and an indicator on whether the product has a price reduction of more than 5%.
The goal of our study is to predict the future sales volumes of each product from each store given each promotion strategy based on historical sales data. Through this prediction procedure, we are able to estimate future purchases, evaluate the influence of promotion strategy for product sales, and potentially recommend the most profitable products to store managers, so the company can make wiser decisions on marketing strategies and inventory planning. For considering the trend of product sales, we aggregate the weekly data into monthly data according to the record time information so that the data contain more than 79.2 million sales records for 132 months from the beginning of 2001 to the end of 2011.
For the proposed method, we classify stores, products, observed time points and promotion strategies into subgroups based on their markets, product categories, month of the year and whether a price reduction is applied, respectively.   Table 2 shows the summary statistics of the data. According to the proposed method, the data can be reframed into monthly third-order tensors by store, product and promotion. Figure 5 provides an illustration of the reframed sale records to clarify the proposed method.
According to the given structure of monthly third-order tensors, the total number of sale  Moreover, we show average sales volumes for three arbitrary categories of products over 132 months in Figure 7, and the prediction interval estimated by the proposed method in Figure 8, while the corresponding figures for the rest of the product categories can be found in the Supplementary Material. In Figure 7, the black lines indicate the observed values, and the red lines indicate the estimated values. The time intervals from the 121st month to the 132nd month show the forecasting performance. The DTRSar has performance similar to the DTRSin, and is therefore not provided here. We notice that the DTRSin can estimate forecasting more accurately. Although the REM can estimate sales volumes sufficiently well on the training set, the forecasting on the testing data has relatively larger biases, while the BPTF has poor performance on both the training and testing sets. Figure   8 provides pointwise prediction intervals for average sale volumes estimators under 95% and 50% nominal coverage probabilities.

Discussion
In this article, we propose a new dynamic tensor recommender system which incorporates time information through a tensor-valued function. A unique contribution of the proposed method is that it can effectively forecast future recommendations at irregular time points.
Technically, the proposed method builds a time-value tensor decomposition model and bor-rows group information from existing time points of the same group for higher forecasting accuracy. Moreover, the proposed method utilizes the polynomial spline method and the weighted least squared method to incorporate time-dependency and intra-cluster correlation into the DRS. In addition, the proposed method is able to provide pointwise prediction intervals based on the established asymptotic property, while existing recommender systems are not equipped with prediction intervals. In theory, we demonstrate that the proposed decomposition achieves asymptotic consistency on prediction and the spline coefficient estimators have asymptotic normality. The proposed method shows numerical advantages compared to existing methods. In real example analysis for the IRI marketing data, the proposed method achieves better performance on forecasting than competitive approaches.

Supplementary Material
The online Supplementary Material contains the other results of real data analysis and all technical conditions and proofs.