With the recent advancement of data collection technologies, large-scale event time data recorded with high temporal resolutions have become increasingly available through various trading, social media, and retail platforms. Analyzing such data is important in understanding user behaviors. One major challenge with such data is that they often have large sample sizes and complex structures. As an example, we analyze a dataset that draws from stock trading transactions (recorded at the second level) by more than 300,000 Chinese trading accounts over approximately three years. Figure 1 shows the estimated overall trading intensity within a transaction day and the ratios of the estimated trading intensities at the account and day levels relative to the overall trading intensity, based on 1,000 sample accounts over 672 days. As shown in the figure, the trading intensities vary considerably within each day, between accounts, and across days. Our main interest is to develop a methodology to investigate the impacts of the di↵erent sources of variations for event time data that are repeatedly observed over time and from multiple accounts, as is the case in our example data. To that end, we view the event times from each account on a given day as a realization from a so-called log-Gaussian Cox process (LGCP, Møller et al., 1998;Simpson et al., 2016). A log-Gaussian Cox process is a Poisson process that has a latent intensity function the logarithm of which is a Gaussian process. We decompose the log-latent intensity as the sum of three independent Gaussian processes, reflecting variations at the account, day, and account-day levels, respectively. To maintain modeling flexibility, we approximate these latent Gaussian processes nonparametrically using the Karhunen-Loève expansion (Watanabe, 1965), in which each process is characterized as a linear combination of orthogonal eigenfunctions (as illustrated in ( 6)). This approach does not require parametric assumptions such as specific forms of the mean and the covariance function of the latent Gaussian process that are commonly used in existing literature (see Yu et al., 2009;Simpson et al., 2016, etc.). Instead, we impose some less restrictive moment and smoothness conditions on the latent Gaussian processes, which are outlined in detail in Section 4. In this sense, the proposed approach o↵ers greater flexibility in modeling complex event time data.
Although the use of the Karhunen-Loève expansion (Watanabe, 1965) in the current setting may appear straightforward, there are nontrivial challenges. In particular, one challenge is how to nonparametrically estimate the covariance function for each latent Gaussian process based on which the eigenfunctions and associated principal component scores can be estimated. Instead of trying to estimate the covariance functions directly, we first propose four nonparametric estimators to estimate di↵erent aspects of the covariance structure, and then define the estimator for each covariance function as a ratio involving a subset of these four estimators. To predict the principal component scores, we develop a conditional likelihood estimation approach. Both procedures are novel compared to their counterparts in the existing literature (see Section 2 for details). Furthermore, we allow correlation among the random principal component scores over time. In contrast, most existing literature assumes the principal component scores to be independent. Although allowing this correlation creates complexity in studying the theoretical properties of the proposed estimators, this complexity is accommodated in our framework (see Section 4 for details).
Our proposed modeling framework is extremely flexible and can be readily extended to model event time data with even more complex structures. In our example data, each event is either a buy or sell transaction. Therefore, we model the observed event times from each account on a given day as a realization from a bivariate point process instead of a univariate point process. To understand the correlation between the buying and selling activities, we extend our proposed modeling framework to a bivariate log-Gaussian Cox process model (see Section 3 for details). In addition to the transaction types, the data contain information on the brokerage branch to which each account belongs. A brokerage branch is a geographic o ce at which an individual establishes a stock trading account, and its location is typically close to the account holder's home address. It may be interesting to study whether the trading activities from accounts in the same branch are similar, possibly due to branch-related factors such as local income levels, economic activity, news, social interaction, etc. To answer this question, we can introduce an additional Gaussian process when decomposing the log-latent intensity function to reflect the branch level e↵ect. In Appendix C, we illustrate how our proposed modeling framework can be easily generalized to this new setting.
There has been a growing interest in using point process models to analyze event time data in the machine learning community (see Zhao et al., 2015;Karimi et al., 2016;Farajtabar et al., 2016Farajtabar et al., , 2017;;Hosseini et al., 2017;Xiao et al., 2017;Xu et al., 2017;Upadhyay et al., 2018;Zarezade et al., 2018, etc.). Most existing work is built upon parametric point process models such as the Hawkes process (Hawkes and Oakes, 1974). Several multivariate Hawkes process models have been proposed to jointly model sequences of observed events in continuous time (see Farajtabar et al., 2016Farajtabar et al., , 2017;;Zarezade et al., 2018;Achab et al., 2018, etc.). In principle, for our example data, the multivariate Hawkes process can be used to directly model potential interactions between trading activity among di↵erent accounts. However, a direct application of the multivariate Hawkes process to our example data is di cult at least for two reasons. First, the multivariate Hawkes process depends on a set of background rate functions that characterizes the overall rate of events for di↵erent accounts and days, which can be rather complex and heterogeneous, as suggested by Figure 1. It can be challenging to estimate such a large number of background rate functions simultaneously. Second, it is unclear how one should define the triggering kernel across di↵erent accounts and days. The triggering kernel is typically a function of the distance between a given time point and a past event time. Because the stock market does not operate on a continuous-time schedule, it is not straightforward to define the distance between two time points from di↵erent days and possibly also from di↵erent weeks.
From a di↵erent perspective, our example data are replicated time series of point processes, in that the trading times of each account form a series of point processes and the di↵erent accounts constitute the replicates. Such data are di↵erent from that typically considered for Hawkes processes, and consequently would require a di↵erent modeling approach. Nevertheless, for the stock trading data under consideration, one can still fit a univariate Hawkes process model to the combined temporal point process from each account, (see Section 2.2 for more details). Our simulation study and the real data analysis demonstrate that such a univariate Hawkes process may fail to capture important aspects of the observed multi-level temporal point process data.
In contrast to the multivariate Hawkes process, the proposed method aims at modeling the correlations between trading activities from di↵erent accounts by identifying the contributing latent factors. For example, under our model (1), trading activities from di↵erent accounts on the same transaction day are correlated due to a common underlying day e↵ect. The latent intensities used in our proposed modeling framework are random functions. Our approach, therefore, connects the literature on log-Gaussian Cox processes and multi-level functional data analysis, where the term "multi-level" refers to the decomposition of the log-latent intensity functions at the account, day, and account-day levels. Even though modeling repeatedly observed random functions is well studied in the functional data analysis literature (e.g., Yao et al., 2005;Hall et al., 2006;Di et al., 2009;Yu et al., 2009;Li et al., 2010Li et al., , 2013;;Chen and Müller, 2012;Chen et al., 2017), existing techniques in classical functional data analysis do not apply to our setting because the latent intensity functions are not observed.
A few recent works consider functional data analysis for point processes (e.g., Bouzas et al., 2006;Illian et al., 2006;Wu et al., 2013;Gervini, 2016;Xu et al., 2020). The modeling techniques used in these work, which ranges from the functional principal component analysis of stochastic density processes (Wu et al., 2013) to the independent component analysis of latent intensity functions (Gervini, 2016), require independent replicates of the event times. However, as seen from Figure 1, the event times from the same account or the same day may be correlated due to shared common trading patterns at the account or day levels. It is unclear whether the existing approaches can be generalized to this more complex setting. In contrast, by using multi-level log-Gaussian Cox processes, we can decompose the log-latent intensity functions of the point processes to meaningfully account for di↵erent sources of variations and to study these variations separately.
Lastly, we note that although there are many work on multivariate point processes (e.g., Van Lieshout and Baddeley, 1999;Baddeley et al., 2000;Lemonnier and Vayatis, 2014;Jalilian et al., 2015;Bacry and Muzy, 2016;Farajtabar et al., 2016Farajtabar et al., , 2017;;Achab et al., 2018), the work is sparse in functional settings with multiple levels. Our proposal, therefore, fills an existing gap in the literature on multivariate point processes.
The remainder of the article is organized as follows. Section 2 introduces the proposed procedure for univariate log-Gaussian Cox processes, which includes the model formulation, estimation, and principal component score prediction. Section 3 extends the procedure to bivariate log-Gaussian Cox processes. Section 4 studies the theoretical properties of the proposed estimators. Specifically, we allow autocorrelations within the random principal component scores in our theoretical framework. Section 5 demonstrates the e cacy of the proposed methods through a simulation study. Section 6 applies the proposed method to the stock trading dataset. All proofs are contained in the Supplementary Material. Throughout the rest of the article, we refer to our proposed approach as the multi-level functional principal component analysis (MFPCA).
To illustrate our model formulation, we describe our model specifications in the context of the stock trading data introduced in Section 1. However, the proposed framework is general and can be used for other applications. One such example is datasets that collect time-stamps of user/consumer activities (e.g., liking or sharing a social media post, making a purchase, posting a review of a service product, etc.) for a large group of users over a certain period. Moreover, a post can be either an original post created by the user or a repost of another user's post, just like each trading event can be either a buying or selling event.
Our proposed method can be applied to investigate the sources of variations in the activity of social media users, and such analyses can provide useful insights into user/consumer activity patterns, which in turn can be used as features in downstream business strategy development such as advertisement placement and content recommendations.
Suppose that we observe daily trading times from n accounts during m days. The trading times observed from account i on day j can be written as N ij = {T ij,l , l = 1, . . . , # {N ij }}, where 0 T ij,1 < T ij,2 < • • • < T ij,#{N ij } T for some given T > 0, and # {N ij } is a nonnegative integer random variable, i = 1, . . . , n, j = 1, . . . , m. We assume that each N ij is generated from an inhomogeneous Poisson process conditional on a latent intensity function ij (t), t 2 [0, T ], where
(1)
In (1), 0 (t) is a deterministic baseline function representing the average daily trading pattern, and X i (t), Y j (t) and Z ij (t) are Gaussian processes characterizing deviations from the baseline for the ith account, the jth day, and the residual deviation, respectively. Model (1) assumes that stock trading activities are a↵ected by both account and day related factors.
The latent process Z ij (t) reflects that di↵erent individuals may adjust their trading behaviors di↵erently in reaction to market conditions or news on a given day. We assume that at any given time t, X i (t), Y j (t) and Z ij (t) are independent normal random variables with mean 0 and bounded variances 2 X (t), 2 Y (t) and 2 Z (t). We further assume that X i (t) and X i 0 (t) are independent for i 6 = i 0 , and that Z ij and Z i 0 j 0 are independent for j 6 = j 0 . However, Y j (t) and Y j 0 (t) can be correlated even if j 6 = j 0 ; this enables us to account for potential temporal dependences over consecutive transaction days.
The motivation for model (1) can be better understood in the context of our example stock trading data. An investor's trading activity on a given day is jointly impacted by various factors which include personal characteristics (e.g., social-economic status), and the overall stock market condition that is a↵ected by the day-related events such as macroeconomic conditions, market-wide news, government policy interventions, etc. An important goal of our modeling approach is to extract account-related contributions to trading activity while controlling for the impacts of day-related factors. Such extraction can help a researcher investigate potential di↵erences in trader behavior. Using model (1), we can accomplish this goal by studying the di↵erences among fitted X i (•)'s for i = 1, • • • , n. On the other hand, exploiting the variations at the day level could help a researcher understand how evolving macroeconomic conditions, government interventions, etc., a↵ect an investor's trading behavior. This can be accomplished by studying how Y j (•)'s evolve over time as j increases from 1 to m. Section 6 contains more detail on these issues.
To the best of our knowledge, there has been no existing method that directly models multilevel temporal point process data as described in Section 2.1. For the stock trading data, one simple strategy is to create a single-level process for each account i by combining N ij , j = 1, • • • , m, sequentially, as illustrated in Figure 2. As a result, the combined process, denoted as N c i , becomes a temporal point process defined on the extended time interval [0, mT ], and can be modeled using existing methods. The first model we consider is the Hawkes process. Denote N c i (B) as the number of events in B ⇢ R from account i, and define the corresponding intensity function as
where H i,t is the event time history up to time t for the account i, i = 1, • • • , n. The intensity function of a Hawkes process takes the form
where µ c i (•) > 0 is the background intensity, and ! i (•) is the transfer function, for i = 1, • • • , n. The transfer functions are typically assumed to be nonnegative, suggesting that past events increase the chance of occurrence for future events, commonly known as the "self-exciting" property. For the stock trading data, it is reasonable to assume that the background intensities are the same for a given account at di↵erent days, which can be approximated by cyclic cubic B-splines as follows:
where B k (•)'s are cyclic B-spline basis functions defined on [0, 1], and btc is the largest integer less than or equal to t. For the transfer function, we use the popular exponential triggering kernel
All parameters can be estimated using standard maximum likelihood estimation (Rizoiu et al., 2017).
The second benchmark model we consider is the Log-Gaussian Cox process with a parametric covariance function (P-LGCP). Specifically, we assume that for each N c i , there exists a latent Gaussian process X c i (t) with a mean 0, and an exponential covariance function
where the background intensity c i (•)'s are also approximated by cyclic B-spline basis functions as suggested in (3). All parameters in model ( 4) can be obtained by the two-step estimation procedure proposed in Waagepetersen and Guan (2009).
A closer look at Figure 2 suggests some potential drawbacks of models ( 2) and (4). For each N c i , the latent process X i (•) appears repeatedly for m days, and hence its impacts can be incorporated to the nonparametric periodic background intensity function (3) for both models (2) and (4). Therefore, variations at the account level can still be captured to a certain degree by modeling N c i 's separately using the benchmark models. However, such a strategy is not suitable for evaluating the impacts of Y j (•)'s. Specifically, for the Hawkes process model (2), the self-exciting property may not completely characterize the clustering patterns because the event occurrences are also triggered by day-related external factors Y j (•)'s, which impact all accounts at any given day j. The same issue exists with the P-LGCP model (4). In our simulation study and real data analysis, both benchmark models fail to capture variations arise from Y j (•)'s. Our simulation study further suggests that they also yield worse predictions of intensities ij (•)'s when the observed point process data are generated from the multi-level model (1). See Sections 5.2 and 6.4 for more details.
In this subsection, we describe the proposed estimation framework for model (1). Following the formula of the moment generating function of a Gaussian random variable, the marginal intensity function can be calculated as
where the expectation is taken jointly over X i (t), Y j (t) and
be the covariance function of X i (t). Suppose that R X (s, t) has the following spectral decomposition: andY k 's and Z k 's are the corresponding orthogonal eigenfunctions. Using the Karhunen-Loève expansion (Watanabe, 1965), we can express the random functions X i (t), Y j (t) and Z ij (t), i = 1, . . . , n, j = 1, . . . , m, as
where ⇠ X ik , ⇠ Y jk , and ⇠ Z ijk are normal random variables with mean 0 and variances ⌘ X k , ⌘ Y k and ⌘ Z k . The expressions in (6) have infinite parameter space dimensionality and are infeasible for estimations in practice. Instead, a common approach is to approximate (6) by using the leading principal components, i.e.,
where p X , p Y and p Z are su ciently large to include dominant modes of variation, and the notation '⇡' indicates potential truncation error when using only finite p X , p Y and p Z . We write
. . , n, j = 1, . . . , m, and refer to them as the functional principal component scores.
In our model, we assume that the functional principal component scores at the account, day and account-day levels are independent of each other, i.e., ⇠ X i , ⇠ Y j and ⇠ Z ij are uncorrelated. We assume that ⇠ X ik and ⇠ X i 0 k 0 are independent if (i, k) 6 = (i 0 , k 0 ), that ⇠ Y jk and ⇠ Y j 0 k 0 are independent if k 6 = k 0 , and that ⇠ Z ijk and ⇠ Z i 0 j 0 k 0 are independent if (i, j, k) 6 = (i 0 , j 0 , k 0 ). However, ⇠ Y jk and ⇠ Y j 0 k can be correlated, i.e., the principal component scores at the day level can be autocorrelated. For example, ⇠ Y jk 's may follow an autoregressive correlation structure, i.e., ⇠ Y jk = P q Y l=1 ✓ lk ⇠ Y (j l)k + " jk , for some positive integer q Y , where " jk 's are independent random errors.
To estimate the eigenvalues and eigenfunctions (⌘
Existing methods (see, e.g. Di et al., 2009) cannot be used to estimate these functions since the intensity functions are unobserved (i.e., latent functions).
The marginal second-order intensity ⇢ ij,i 0 j 0 ,2 (s, t) = E[ ij (s) i 0 j 0 (t)] can be shown as
for i, i 0 = 1, . . . , n, j, j 0 = 1, . . . , m. To ease presentation, we first focus on the case where Y j (s) and Y j 0 (t) are independent if j 6 = j 0 . We point out that our method is valid even when this assumption does not hold. As a matter of fact, our theoretical investigation is carried out without this independence assumption. Under this assumption, it follows from (5) and the moment generating function of a Gaussian random variable that s, t)] and D(s, t) ⌘ ⇢(s)⇢(t). Note that A(s, t), B(s, t) and C(s, t) contain information about the correlation at the account-day, account and day levels, respectively. The covariance functions R X (s, t), R Y (s, t) and R Z (s, t) can be easily calculated once we obtain A(s, t), B(s, t), C(s, t) and D(s, t). Next, we develop estimators for these functions, which are later used to derive estimators for the covariance functions.
We first consider the estimation of ⇢ ij,i 0 j 0 ,2 (s, t). By Campbell's theorem (e.g., Møller and Waagepetersen, 2003), it holds for any measurable function
where the expectation is over the point processes N ij and N i 0 j 0 . Let K(•) be a kernel function and K h (t) = h 1 K(t/h). We can estimate ⇢ ij,i 0 j 0 ,2 (s, t) under the four di↵erent scenarios in
dx is an edge correction term. To understand the construction of these estimators, consider b A h (s, t) as an example. By defining
If the bivariate function A(•, •) is smooth and h is su ciently small such that A(u, v) ⇡ A(s, t) for any (u, v) in the small neighborhood around (s, t) defined by the function
for s, t 2 [0, T ]. As suggested by Theorem 1 in Section 4.1, one may need to use di↵erent bandwidths h x , h y and h z for estimating R X (s, t), R Y (s, t), R Z (s, t) when n, m are of di↵erent magnitudes.
From the estimated covariance functions, the eigenvalues and eigenfunction can be estimated using methods proposed in Yao et al. (2005)
where b X k 's are subject to constraints
)'s can be estimated similarly. The validity of solving (11) for eigenvalues and eigenfunctions is ensured by the uniform convergence of the covariance function b R W,h (s, t) to R W (s, t) for W = X, Y, or Z under mild conditions, as given by Theorem 1 in Section 4.1. We refer the reader to the proof of Theorem 2 in Yao et al. (2005) for more technical details.
Remark 1 The computational cost of the estimators given in (10) is of order O(nm), or more precisely, of the order of the total number of events in all N ij 's. We use b A h (s, t) as an example to illustrate this point. Note that
.
The computational cost to calculate each of the three single summations in the above is of the order of the number of events in N ij . Thus, the computational cost to obtain b A h (s, t) is of the order of the total number of events in all N ij 's, which is of order O(nm). Similarly, we can show that the same is true for the other three terms. Therefore, the overall computational complexity of (10) is of order O(nm), which is confirmed by our simulation studies in Section 5.1. See Figure 5 for more details.
Having estimated the eigenvalues and eigenfunctions, we next predict the functional principal component scores. Markov Chain Monte Carlo based approaches have been proposed for predicting principal component scores (Di et al., 2009). However, such approaches are not applicable in our setting because our model cannot be formulated as a linear mixed model (see Di et al., 2009). In this section, we develop a new procedure based on conditional likelihood. We will focus on predicting the principal component scores at the account and day levels, i.e., ⇠ X i 's and ⇠ Y j 's. We give an approach to estimate ⇠ Z ij 's in Appendix B. To predict the functional principal component scores for X i (t), define
In the stock trading data, N X i and N X i are the aggregated trading times over m days for account i and for all accounts excluding account i, respectively. Conditional on the latent function X i (•), the marginal intensity functions of N X i and
where ij (t) is as defined in (1), the first expectation is taken over Y j (t)'s and Z ij (t)'s, and the second expectation is taken further over X i 0 (t)'s with i 0 6 = i. Using the approximation of X i (t) in (7) with estimated eigenfunctions b X k (t)'s and conditional on the presence of an event at t, the probability for this event to be from account i is
where b X k (t)'s are obtained from solving (11) and b 2 X (t) = RX,h (t, t) is obtained from (10). We then predict ⇠ X i by maximizing the following conditional likelihood:
Similarly, the functional principal component scores for Y j (t) can be obtained by maximizing the conditional likelihood
where N Y j and N Y j are the aggregated trading times over n accounts on day j and for all days excluding day j, respectively, and
.
In this section, we propose a goodness-of-fit measure for model (1). To demonstrate the main idea, we consider the aggregated trading times over m days for di↵erent accounts, denoted as N X 1 , • • • , N X n , as defined in Section 2.5. Let 0 = r 0 < r 1 < • • • < r L = T be a sequence of time points that form a partition of the interval [0, T ], and denote by #{N X i }(r l 1 , r l ) the total number of events from N X i that fall between the sub-interval (r l 1 , r l ], l = 1, • • • , L. Focus on the sub-time interval (r l 1 , r l ], the proportion of observed events in [ n i=1 N X i (r l 1 , r l ) that are from the ith account, i.e. from N X i (r l 1 , r l ), can be estimated as
When the number of days m is large, pi (r l 1 , r l ) under model ( 1) is a good estimator for
where ⇢X i [t|X i (•)]'s are as defined in (12). Note that ⇢X
⇤ due to (5) and ( 12), i = 1, . . . , n, where ⇢(t) is defined in (5). Using the estimated model (1) with a bandwidth h, we can then approximate the p i (r l 1 , r l )'s as
Note that we have deliberately emphasized the dependence of pi (r l 1 , r l ; h) on the bandwidth h. The goodness-of-fit of model ( 1) at the X level within the sub-time interval (r l 1 , r l ] can be measured by the discrepancy between the two discrete probability distributions b
We propose the following goodness-of-fit measure based on the Kullback-Leibler divergence
for l = 1, • • • , L. In practice, the function D X (r; h) is a step function on [0, T ] that can be used as a graphical tool to check whether model (1) with the bandwidth h fits the data adequately well at the X-level. A goodness-of-fit measure at the Y level, denoted as D Y (r; h), can be defined in the same way using the aggregated data over n accounts. Finally, our numerical studies suggest that the choice of partition 0 = r 0 < r 1 < • • • < r L = T is of limited importance and we recommend using equal quantiles of the pooled event times
[ n i=1 N X i with L = 20. Such a choice ensures that there are 5% observed event times within each sub-interval (r l 1 , r l ) for l = 1, • • • , L.
Next, we describe a procedure to select the bandwidth in (9) based on the goodness-of-fit measure proposed in Section 2.6. As suggested by Theorem 1, the optimal bandwidths at X, Y, Z levels might be di↵erent, and they depend on n and m. As such, our proposed method chooses di↵erent bandwidths for R X (•, •), R Y (•, •) and R Z (•, •), respectively. We first consider the bandwidth selection for R X (•, •). Randomly partition the
The main idea is to use data in S k = [ k 0 6 =k S k 0 to fit model (1) and then measure its goodness-of-fit to the hold out data in S k using the approach proposed in Section 2.6. The cross-validation score can then be defined as
where b P X,S k l and PX,S k l,h are as defined in (18) using data included in S k and S k , receptively.
The optimal bandwidth ĥx is then chosen by minimizing CV X (h). Following the same procedure, we can define the cross-validation score CV Y (h) using data aggregated over all accounts, N Y 1 , • • • , N Y m , as defined in Section 2.5 and choose the optimal ĥy accordingly. The bandwidth selection for R Z (•, •) cannot be carried out in the same way since we can not predict ⇠ Z ij 's consistently. Alternatively, we simply choose the ĥz for R Z (•, •) as the one that provides good fits at both X and Y levels by minimizing the added cross-validation score
Our simulation studies suggest that such a choice yields satisfactory performances.
In the stock trading data, we observe two types of events from each account, i.e., buy or sell trades. Therefore, we can view the observed event times as realizations from bivariate point processes. By modeling the underlying bivariate point processes, we may gain insight into the temporal covariation of accounts' buying and selling activities.
Let N 1,ij and N 2,ij denote the buying and selling times observed from account i on day j, i = 1, . . . , n and j = 1, . . . , m. As in Section 2, we assume that N r,ij is generated from an inhomogeneous Poisson process conditional on a latent intensity function r,ij (t), where
Similar to the univariate case, r,0 (t) is a fixed baseline function representing the average daily pattern, X r,i (t) and Y r,j (t) are random functions reflecting deviations from the baseline at the account and day levels, respectively, and Z r,ij (t) is the residual deviation, r = 1, 2.
We assume that at any given t, X r,i (t), Y r,j (t) and Z r,ij (t) are independent Gaussian random variables with mean 0 and variances 2 r,X (t), 2 r,Y (t) and 2 r,Z (t), r = 1, 2. Similar to the development in the univariate case, we can write X r,i (t), Y r,j (t) and Z r,ij (t) using the Karhunen-Loève expansion as:
where p (r)
Y and p (r)
Z are su ciently large positive integers, ⇠ X r,ik , ⇠ Y r,jk and ⇠ Z r,ijk are normal random variables with mean 0 and variances ⌘ X r,k , ⌘ Y r,k and ⌘ Z r,k , and X r,k (t), Y r,k (t) and Z r,k (t) are orthogonal eigenfunctions. Define ⇠ X r,i = (⇠ X r,i1 , . . . , ⇠ X r,ip
) T , r = 1, 2. At each r, we assume the same correlation structure within and between ⇠ X r,i , ⇠ Y r,j and ⇠ Z r,ij as in the univariate case. To introduce dependence between two point processes, we further assume that (⇠
) each follow a multivariate Gaussian distribution with covariance matrices ⌃ X , ⌃ Y and ⌃ Z , respectively. Under such assumptions, we can write
where
), W = X, Y, Z and r = 1, 2.
The eigenvalues and eigenfunctions (⌘ X r,k , X r,k (t))'s, (⌘ Y r,k , Y r,k (t))'s and (⌘ Z r,k , Z r,k (t))'s can be estimated using the approach detailed in Section 2.4; the principal component scores can be predicted using the approach described in Section 2.5. The primary goal in this section is the estimation of the covariance matrices ⌃ X 12 , ⌃ Y 12 and ⌃ Z 12 .
We first describe cross covariances estimation at the X level. To that end, first define
This implies that
Because X 1,k (s) and X 2,k 0 (t) can both be estimated, an estimate of X 12,kk 0 can be obtained if Q X (s, t) can be estimated. This can be achieved by generalizing the approach used to estimate the covariance function R X (s, t) in Section 2.4. The cross covariances at the Y and the Z levels can be estimated similarly. The detailed formulas are in Appendix A.
In this section, we investigate the theoretical properties of our proposed estimators when both the number of user accounts n and the number of days m increase. We show that the covariance and cross-covariance function estimators suggested in ( 10) and ( 26) converge uniformly in probability to their theoretical counterparts as n, m ! 1. Consequently, all eigenfunctions and eigenvalues at di↵erent levels as well as correlations among functional principal component scores of buying and selling activities can be consistently estimated. We show that with su ciently large m and n, the functional principal component scores for each account, i.e., ⇠ X i and each day, i.e., ⇠ Y j , can be consistently estimated. Asymptotic convergence rates are derived for all consistent estimators. All proofs are collected in the Supplementary Material.
We start our investigation on the asymptotic properties of the covariance function estimators in (10). Before we state our theorem, we first introduce some regularity conditions.
[C1] There exist constants c 0 > 0 and
C 0 for any s, t 2 [0, T ] and i, i 0 = 1, . . . , n; j, j 0 = 1, . . . , m.
[C2] There exists a constant
[C3] Assume that (a) h x ! 0, nh 4
x ! 1 and nh 6
x < 1; (b) h y ! 0, mh 4 y ! 1 and mh 6 y < 1; and (c)
Assume that there exists a constant C M such that, with a probability tending to 1,
for any s 1 , s 2 , t 1 , t 2 2 [0, T ], and W = X or Y .
Condition C1 imposes some constraints on the random daily intensity functions ij (•)'s. Specifically, condition C1(a) requires that, on average, there should be only a finite number of transactions for any user on a given day. Condition C1(b) is satisfied if the baseline intensity 0 (•) is bounded from above and below since we have assumed that sup s2[0,T ] 2 W (s) < 1 for W = X, Y, and Z. From a practical point of view, this means that the average user trading frequency in a stock market should be finite and bounded away from 0 at any time of day. Condition C1(c) imposes some mild smoothness conditions on 0 (t) and the covariance functions R W (s, t), W = X, Y , and Z, so that they have bounded absolute partial derivatives. Condition C2 essentially assumes that there are only short term longitudinal correlations among trading intensities in the stock market. Conditions C3-C4 are technical conditions on the choice of bandwidth and kernel functions that are analogous to those in, for example, Yao et al. (2005). Condition C5 imposes some mild restrictions on the smoothness of the sample paths of the latent stochastic processes X i (•)'s and Y j (•)'s, and they closely resembles the commonly employed stochastic equicontinuity assumption in the literature (Newey, 1991) The following theorem states that with the above conditions, the covariance estimators in (10) are uniformly convergent in probability to the true covariance functions.
Theorem 1 Under conditions C1-C5, we have that as n, m ! 1,
The proof is given in the Supplementary Material. Theorem 1 states that the convergence rates of b R X,hx (s, t) and b R Y,hy (s, t) are determined by the number of replicates at the X level (i.e., n) and the Y level (i.e., m), respectively. However, the convergence rate of b R Z,hz (s, t) is determined by n and m jointly. These observations are confirmed by our simulation studies in Section 5.1.
The proof of Theorem 1 utilizes similar technical devices used in Yao et al. (2005), but several new challenges need to be addressed. Firstly, unlike Yao et al. (2005), the intensity functions are latent functions that are not directly observable. Therefore, we need to study properties of a series of carefully designed random sums given in (9), instead of simply smoothing observed functional data. As a result, new conditions on the observed point processes need to be imposed as suggested in condition C1. Secondly, Theorem 1 investigates multi-level latent functions while Yao et al. (2005) only studies single-level functional data. To the best of our knowledge, there are limited rigorous theoretical investigations on the asymptotic properties of multi-level functional data analysis despite its popularity. Therefore, our work makes a useful contribution to the literature in this respect. Finally, the functional data studied in Yao et al. (2005) are independent random trajectories while the observed point processes in our work are not. For example, all N ij 's for di↵erent i's are dependent due to the common component Y j in their latent intensities defined in (1). In addition, we also allow dependence among the latent processes Y j (t)'s in Theorem 1 as specified in condition C2. The involvement of dependence further complicates our analysis.
The following corollary is a direct consequence of Theorem 1.
Corollary 1 Under conditions C1-C5, we have that, as n, m ! 1,
The proof of Corollary 1 follows readily from the proof of Theorem 2 in Yao et al. (2005) and is thus omitted.
Theorem 1 and Corollary 1 ensure the accuracy of the estimated covariance functions and their corresponding eigenfunctions. It is therefore reasonable to expect the conditional likelihoods ( 13) and ( 14) to yield consistent estimators ⇠ X i and ⇠ Y j , which characterize the trading pattern at the account and day levels, respectively. Denote k • k as the Euclidean norm. Our theoretical investigation relies on the following conditions.
[S1] E
[S5] The number of functional principal components, p W , satisfies (a) (n 1 +m 1 )log(m+
Conditions S1-S5 are technical conditions on the latent Gaussian processes. Condition S1 imposes mild moment conditions on the latent Gaussian processes. Condition S2 is trivially true if all eigenfunctions in X (•) and Y (•) are uniformly bounded over the domain [0, 1]. Condition S3 ensures that the Hessian matrix of the conditional likelihoods ( 13) and ( 14) are strictly positive definite with a probability approaching one. Condition S4 is only needed for Y j (•) so that there are only a "small" number of Y j 0 (•) strongly correlated with Y j (•), analogous to condition C2. Condition S5 requires that the decaying rates of eigenvalues of covariance functions should be su ciently fast as p W grows and that p W should not be too large compared to sample sizes n, m. These additional conditions are needed as latent random functions are not observed in our framework. Moreover, all components in ⇠ X i or ⇠ Y j need to be estimated simultaneously. This is di↵erent from classical functional data analysis where the principal component scores can be estimated separately, see, e.g., Yao et al. (2005); Chen and Müller (2012).
Theorem 2 Let b ⇠ X i and b ⇠ Y j be maximizers of the conditional likelihoods (13) and (14), respectively. Under conditions C1-C5 and S1-S5, and assuming that all eigenvalues of covariance functions R X (s, t) and R Y (s, t) are of multiplicity 1, we have that as p X , p Y , n, m ! 1,
The proof is given in the Supplementary Material.
Theorem 2 suggests that the convergence rates of b ⇠ X i and b ⇠ Y j are controlled by (i) the estimation errors of the covariance functions and their eigenfunctions; and (ii) the truncation errors of the Karhunen-Loève expansion in (7), which are reflected in terms e p X and e p Y . Note that the lengths of the vectors ⇠ X i and ⇠ Y i increase as p X and p Y grow, and hence the overall convergence rate in Euclidean norm k • k will be slower by a factor of p
To study the convergence rate of the cross-covariance function estimators given in ( 26), we need some additional assumptions on the marginal cross second-order intensity functions and the cross covariances among the day level latent processes.
[C1*] There exist constants c 0 > and
Assume that there exist a constant C M such that, with a probability tending to 1,
Theorem 3 Under conditions C1*-C3* and C3-C4, we have that, as n, m ! 1,
The proof of Theorem 3 is almost identical to that of Theorem 1 and thus is omitted. The only di↵erence is that one has to impose three additional assumptions C1 ⇤ -C3 ⇤ to regulate the strengths of dependence between the two types of point processes and smoothness of latent processes.
The following corollary is an immediate consequence of Theorem 3.
Corollary 2 Under conditions C1*-C3* and C3-C4, if we further assume that all eigenvalues of covariance functions R r,W (s, t), W = X, Y, or Z, are of multiplicity 1, we have that as
where the convergence is entry-wise for the matrix.
The proof is a straightforward application of the Slutsky's Theorem combining the results of Corollary 1 and Theorem 3. We thus omit the proof.
In this section, we demonstrate the e↵ectiveness of the proposed approach through simulation studies. For the univariate log-Gaussian process, the data are simulated from the univariate point processes model presented in Section 2.1 using the following intensity model:
and simulate principal component scores ⇠ X ik , ⇠ Y jk and ⇠ Z ijk as independent normal random variables with means 0 and variances equal to the eigenvalues
sin(4⇡t)}, which are used to mimic those in our real data. Throughout this section, the Epanechnikov kernel function is used in all simulation studies and summary statistics based on B = 500 simulation runs are calculated.
We first evaluate the accuracies of the proposed MFPCA procedure in estimating the eigenvalues and eigenfunctions for univariate point processes. In this subsection, we fix
2 and allow n, m to vary. The bandwidths ĥx , ĥy , and ĥz are selected by the data-driven procedure proposed in Section 2.7.
Figure 3 shows the estimation accuracies of the eigenfunctions based the mean absolute deviation (MAD) defined as
where b W k (t) (b) 's are the estimated eigenfunctions in the bth simulation. Figure 3 (a) and Figure 3 (d) show that the estimation accuracies of X k (t)'s improve as n increases, but not necessarily so as m increases. Figure 3 Next, Figure 4 illustrates the estimation accuracies of the eigenvalues using the relative root mean squared errors (RMSE) defined as
The results are similar to our observations in Figure 3 and thus we omit further discussion. 13)-( 14) and the true ones under various settings. It is seen that at both the X and Y levels, the correlation coe cients are generally rather high in all case scenarios, suggesting excellent prediction performances. We also observe that compared to the first PC scores, the second PC scores are generally more di cult to predict. Finally, Figure 5 (c) and (f) illustrate the averaged CPU time (in seconds) for estimating the covariance functions R W (s, t), W = X, Y, Z, with a fixed bandwidth. All simulation runs were carried out in the software R on a cluster of 100 Linux machines with a total of 100 CPU cores, with each core running at approximately 2 GFLOPS. As discussed in the Remark 1 in Section 2.4, the computational costs of the proposed estimation method indeed grows linearly with m and n and is of the order O(nm).
In this subsection, we investigate the benefits of using the two-level model compared to benchmark models suggested in Section 2.2, namely, the Hawkes process model and the parametric Log-Gaussian Cox process model (P-LGCP). In addition, we consider a singlelevel model assuming that the point process N ij has the latent intensity ij (t) of the form
We simulate the data from model ( 23) with a fixed
Under this setting, the two-level model reduces to the single-level model when q = 0, and the di↵erence between these two models becomes larger as q grows. All other settings remain the same as those in Section 5.1. To fit the single-level model, we apply the method proposed in Wu et al. (2013) to the aggregated point processes
, to predict the random intensity function in (24). For comparisons between various fitted models, we consider three metrics: (1) the overall mean absolute deviation: OMAD =
(2) the overall goodness-of-fit at the X level: Fit X = R 1 0 D X (r; h)dr with D X (r; h) defined in (18);
(3) the overall goodness-of-fit at the Y level: Fit Y that is similarly defined as Fit X .
For the Hawkes and the P-LGCP model, the estimated background intensity μi (•)'s as suggested in (3) are used as ˆ ij (•)'s in the definition of OMAD. The summary statistics based on B = 500 simulation runs are summarized in Figure 6. When computing Fit X and Fit Y for the Hawkes and the P-LGCP model, the model-based probabilities defined in ( 16) do not have closed-forms as given in (17). Therefore, for each simulation run, we numerically approximate the model-based probabilities by the average of the empirical probabilities (15) computed using 100 replicate datasets simulated from the fitted models. Figure 6 (a) shows that as q increases, the estimation accuracies of ij 's using the single-level model as well as the two benchmark models become much worse than the two-level model. However, the single-level model can indeed still capture the variations at the X level as indicated by Figure 6 (b), where the averaged Fit X for both single-level and two-level models are well within the 95% quantile bands of Fit X (Truth).
In contrast, while the Hawkes and the P-LGCP model can capture the general trend of the variations at the X level, their performances are significantly worse than those of the other two models. One possible explanation is that the parametric assumptions imposed on the benchmark models may not be suitable for the underlying data, as we have discussed in Section 2.2. Since the single-level model ( 24) completely ignores the variations at the Y level, it essentially assumes pj (r l 1 , r l ; h y ) = 1/m for j = 1, • • • , m in the definition of Fit Y , which provides rather poor goodness-of-fit at the Y level as q grows. The Fit Y 's computed using simulated data from the fitted Hawkes and P-LGCP models have almost the same trend as q increases in Figure 6, suggesting that both models fail to capture variations at the Y level. In summary, when data are generated from the two-level model ( 23) but fitted with a single-level model ( 24) or benchmark models given in Section 2.2, variations at the X level can be captured to a certain degree while variations at the Y level may be poorly characterized. This leads to rater inaccurate estimates of the underlying latent intensity functions ij (t)'s. ) is computed by replacing pi (r l 1 , r l ; h) in ( 18) with the theoretical probability (16).
Next, we evaluate the estimation performance for the bivariate point process model presented in Section 3.1. The same eigenvalues and eigenfunctions are used as in the univariate model ( 23). The principal scores are generated from the covariance structure in ( 22), with
The example data to which we apply our methodology contains stock market buy and sell transactions from 331,599 investor accounts at a leading brokerage house in mainland China during the period covering January 4th, 2007 to September 30th, 2009. These accounts are from 98 branches distributed over 25 provinces in the central and eastern coastal areas of mainland China.
As background on the Chinese stock market, during 2006-2017 its market capitalization grew more than five-fold to become the second-largest in the world (Carpenter and Whitelaw, 2017). Among the world's ten largest stock exchanges, two are operated in mainland China (the Shanghai Stock Exchange and the Shenzhen Stock Exchange)foot_0 . The Chinese market began to boom in November 2006 (Andrade et al., 2013), just before the start of our sample period. The Shanghai Stock Exchange (SSE) Composite Index doubled from around 2,000 points to more than 4,000 points in less than six months. Many observers characterized the market as experiencing a bubble, and in an attempt to cool down the market, on May 30th, 2007, the Chinese government tripled the transaction tax on trades from 0.1% to 0.3%.foot_1 Although the market experienced several major sell-o↵s following the tax increase, the SSE market index recovered and kept increasing until reaching a historic high of 6,124 points on October 16th, 2007. Afterward, the market index plunged and in April 2008 the transaction tax was reduced to its prior level of 0.1%. The market continued to decline, with the SSE index in mid-September 2008 down approximately 70% from its October 2007 high. To stimulate buying, on September 18th, 2008, the Chinese Government announced the abolition of the tax on buying (maintaining it only for sales) among other measures to bolster the market. The market gradually recovered and the SSE composite index closed at 2,779 points on September 30th, 2009.
The period we study contains 672 trading days, and on average, each account in our dataset has 150 buy or sell transactions during this period. Each investor account has a randomly-assigned identification number and there is insu cient information for us to discern actual investor identity. There is a large variation in the trading frequencies across accounts where the total number of transactions ranges from 0 to 874,263. In our analysis, we have removed accounts with fewer than 30 or more than 2,000 total transactions, which results in retaining approximately 47% of the accounts in the broader dataset, or equivalently, a total of 157, 203 accounts. We do not consider accounts with fewer than 30 transactions because it would be di cult to estimate the account-level e↵ect for such accounts. Accounts with more than 2,000 trades are likely institutions that submit trades on behalf of other investors and account for only around 0.2% of the total data. Our results are not materially a↵ected by changing the minimum cuto↵ to 10 or by changing the maximum cuto↵ to 1,000.
We next discuss the application of our proposed model to the stock trading data, interpreting what is learned where appropriate. In subsection 6.5, we discuss potential further applications for interested finance researchers.
We first apply our proposed MFPCA approach for univariate point processes to study the variations of an investor's buying and selling activities separately. The intensity functions of the processes underlying the buying and selling activity are modeled using the two-level models (20) with p X = p Y = p Z = 5 and bandwidths for the Epanechnikov kernel at di↵erent levels selected by the data-driven method proposed in Section 2.7. Figure 7 shows the first three eigenfunctions at the account (X), day (Y ), and account-day (Z) levels for both the buying and selling processes. The first three eigenfunctions at the X, Y and Z levels explain about 96.30%, 93.44%, and 67.37% of the variation at the respective levels for the buying process, and 96.62%, 93.25%, 67.44% of the variation for the selling process.
In Figure 7, we can see that the eigenfunctions for the buying and selling processes are very similar. In what follows, we therefore focus on interpreting the eigenfunctions for the buying process alone. Note that the first eigenfunctions of all three levels are flat. This suggests that the overall buying frequency is the most important factor in explaining the variability at each level. At the account level, the second eigenfunction decreases throughout the day except for a small jump at the beginning of the third trading hour. This characterizes trading accounts that display buying activity that steadily increases (or decreases) throughout the day. The third eigenfunction at the account level is approximately flat at zero until an abrupt jump around mid-day; the eigenfunction then decreases from mid-day to the end of the day. This characterizes trading accounts that are more active in the afternoon and whose afternoon trading activities follow an increasing or decreasing pattern. At the day level, the second eigenfunction shows a roughly decreasing trend in the morning. This characterizes days that have increasing or decreasing morning trading activity.
At the account-day level, the second eigenfunction is relatively steady in the morning and afternoon with a sharp drop around noon. This characterizes contrasting morning versus afternoon trading behaviors. The third eigenfunctions at both the day and the account-day levels rise until mid-day and then decrease. This characterizes trading behavior that is concentrated around mid-day.
Figure 8 shows scatter plots of the principal component scores in the first three directions for the buying process at the account level. We can see from these plots that compared to accounts with lower scores in the first direction, those with higher scores tend to have a smaller variation in the second and third directions. Note that a larger principal component score in the first direction indicates higher overall buying frequency. Thus, accounts with higher buying frequency tend to buy more evenly throughout the day than with the specific patterns as described by the second and third directions. Richards and Willows (2018) find that only a small proportion of investors with certain common characteristics trade frequently, which may partially explain the homogeneity of their intra-day trading patterns we observe. However, there has been little research in the finance literature that studies cross-sectional di↵erences in time-of-day trading preferences between investors with higher trading frequencies and those with lower trading frequencies, as suggested by Figure 8. We comment on this finding later below.
Figure 9 shows principal component scores in the first two directions for the buying process at the day level. Because the first two eigenfunctions characterize the overall trading frequency and morning trading activities at the day level (see Figure 7), a larger score in the first direction indicates higher trading frequency, whereas a larger positive (negative) score in the second direction indicates increasing (decreasing) morning trading activity during that day. From the first plot, we can see that the trading frequency increases from the start of our observation window until May 30th, 2007, when the transaction tax is tripled. At that point trading frequency begins to decline, showing that the tax increase was successful in dampening trading activity as the government desired. This finding is also consistent with the findings in transaction tax research (Umlauf, 1993;Mannaro et al., 2008).
Despite noticeable variability, trading frequency continues to decline until September 19th, 2008, the first trading day after the buy-side transaction tax is eliminated. Thereafter, again consistent with the implications of prior research on transaction taxes, trading frequency begins a new trend of increasing intensity. In the second plot, we observe a period of overall decreasing morning trading activity following the tripling of the transaction tax on May 30th, 2007. The largest scores are achieved around September 18th, 2008 when the buy-side transaction tax is abolished, coupled with low scores in the first direction. This implies low trading activity overall that is concentrated early in the trading day. On September 19th, almost all the stocks reach the maximum allowable 10% daily increase shortly after the opening of the market, and there were very few trading activities after that.foot_2 This is consistent with the values of the scores in the first two directions. The scores in the second direction are generally small in the other periods.
Figure 8: Scatter plots of buying scores from di↵erent directions at the account level based on 10000 randomly selected accounts. From left to right: scores of directions 1 vs 2, 1 vs 3 and 2 vs 3.
To study the temporal covariation of an investor's buying and selling activity, we apply our proposed MFPCA approach for bivariate point processes. Table 2 shows the correlation matrices of scores from the buying and selling processes together with the associated standard error estimates obtained by bootstrapping the predicted scores at the account and day levels. The bootstrap procedure is admittedly an ad hoc approach but can nevertheless provide some measure of uncertainty for these estimates. It is not used for the account-day level, however, because of the di culty in estimating the scores reliably at that level (some accounts have too few transactions). Based on the results, we observe strong positive correlations between scores in the same direction and at the same level, i.e., there are large diagonal elements in each of the three correlation matrices. These correlations can be interpreted easily given the interpretations of the associated eigenfunctions. For example, the positive correlation between the first buying and selling directions at the account level (0.9643) suggests that accounts with higher buying frequency also have higher selling frequency. This high correlation indicates that most trades are security selection rather than asset-class-allocation trades, i.e., investors are typically rebalancing an existing stock portfolio rather than increasing or decreasing its overall size. Because it is very unlikely that the typical Chinese individual investor (as opposed to institutional investor) has access to superior information about the stocks they trade, this finding is consistent with the finance literature's documentation that individual investors, to their detriment in terms of performance and trading costs, are overconfident and trade too frequently (Odean, 1999;Barber andOdean, 2000, 2001;Barber et al., 2008). The o↵-diagonal elements in the correlation matrix can also o↵er insights into the correlation between buying and selling activity. For example, at the account level, the negative correlation between the second buying and the third selling directions (-0.1418) implies that accounts with decreasing buying frequencies in a transaction day are also more likely to sell earlier in the day than later. For the day level, the positive correlation between the second buying and the third selling directions (0.3285) implies that transaction days with decreasing buying activity in the morning tend to have higher selling activity around the middle of trading day.
⌃ X ⌃ Y ⌃ Z 0.9643 -0.0005 0.-0.0244 0.7873 0.1053 -0.0167 0.8110 0.1913 -0.0287 (0.1062) (0.0128) (0.0182) (0.0134) (0.0929) (0.0767) 0.0108 0.3982 -0.1418 0.0884 0.4585 0.3285 -0.1636 0.5227 -0.2671 (0.0127) (0.1157) (0.0140) (0.0835) (0.0896) (0.1152) -0.0180 -0.0727 0.7690 -0.0784 -0.293 0.7600 -0.0421 0.2415 0.5936 (0.0130) (0.0335) (0.0560) (0.0425) (0.0928) (0.0234) Table 2: Correlation matrices between scores (⇠ B , ⇠ S ) from the buying and selling processes at each level. Each matrix includes the first three directions. Bootstrapped standard errors are shown in the parentheses.
Finally, we check the goodness-of-fit of the two-level model (1) to the stock trading data.
The results for the buying and selling processes are rather similar, so we only present some diagnostic plots for the buying process. Figure 10 (a) and (c) illustrate the goodness-of-fit measures computed on a random subset of 10, 000 accounts applying the definition (18) at the account (X) and day (Y ) levels with di↵erent numbers of principal components. For both account and day levels, we observe that the model estimated using only one principal component is not able to capture all of the variation across di↵erent accounts or days. As the number of principal components increases to 5, however, the model fits become significantly better for both levels.
As a comparison, we compute the same goodness-of-fit measures for the Hawkes and the P-LGCP model suggested in Section 2.2, with the model-based probabilities numerically evaluated by the average of the empirical probabilities (15) based on 500 replicate datasets simulated from the fitted models. In Figure 10(a), we can see that at the account level, the Hawkes process model performs much worse than all LGCP models. The P-LGCP model performs similar to the MFPCA model with p X = 3 principal components, although its goodness-of-fit at the boundaries is relatively worse. Overall, the MFPCA model with 5 principal components shows uniformly superior performance and achieves rather small KLdivergence, indicating that the model captures account level variation su ciently well. At the day level, Figure 10(d) suggests that both benchmark models perform similarly, yielding model fits that are significantly worse than those of the MFPCA model. This is expected due to the reasons discussed in Section 2.2. The KL-divergence values for the MFPCA model with 5 principal components are also rather small, especially when compared to those at the account level, which indicates an adequate fit at the day level. have found about trading behavior. There are, however, unresolved questions in the finance literature that could benefit from the methodology we propose.
There is a large literature in behavioral finance and social sciences concerning interpretations of trading behaviors. One topic of prior research is how trading activity di↵ers due to investor characteristics such as age, gender, or geographic location (e.g., Odean, 1999;Barber and Odean, 2001;Samanez-Larkin et al., 2010;Henninger et al., 2010, etc.). Another is how information that induces trading disseminates across a large population (e.g., Feng and Seasholes, 2004). However, most existing work is based on aggregated trading activity (e.g., the total number of trades in a day), and thus fails to gain potentially important insights from intraday trading patterns. Furthermore, when studying the impacts of account-related characteristics on trading behavior, trading statistics based on aggregated data across time may fail to account for the confounding e↵ects of day-level factors. As data that include intraday trading information become more available, the proposed MF-PCA framework provides a powerful tool with which to decompose the variation of such trading data in meaningful ways and to summarize high-frequency trading activity using multiple principal component scores. The PC scores provide a multi-dimensional view of intraday trading activity, which in turn allows researchers to more thoroughly investigate and describe the factors that determine trading decisions.
A second line research our proposed method can benefit is that which attempts to exploit abnormal patterns in trading activity to devise successful trading strategies (e.g., Andrade et al., 2008). Standard outlier detection tools can be applied to the accountlevel principal component scores our methodology incorporates to identify investors with abnormal intraday trading activity. This may lead to improvements in trading models by identifying more predictive trading signals. At the day level, it may also be interesting to monitor the daily trading status of the stock market to predict the formation of asset price bubbles (e.g., Andrade et al., 2013), which may better enable policymakers to implement countermeasures. Such monitoring could be accomplished by applying change point detection tools (e.g., Polunchenko and Tartakovsky, 2012) to the multivariate time series of principal component scores at the day level to detect abnormal changes in the stock market, which might indicate bubble formation in its early stages.
In view of ( 28
See the World Federation of Exchanges (https://www.world-exchanges.org).
See Andrade et al. for evidence supporting a stock market bubble during this time-period, as well as further discussion of the tripling in the transaction tax.
During the time period we study, the Chinese markets impose a maximum 10% increase or decrease in a stock's price. Trading in the stock closes for the rest of the day if that limit is reached, which does happen on occasion.
2007-01-01 2007-08-29 2008-04-25 2008-12-21 2009-08-18