Multivariate Gaussian and Student-t process regression for multi-output prediction

Gaussian process model for vector-valued function has been shown to be useful for multi-output prediction. The existing method for this model is to reformulate the matrix-variate Gaussian distribution as a multivariate normal distribution. Although it is effective in many cases, reformulation is not always workable and is difficult to apply to other distributions because not all matrix-variate distributions can be transformed to respective multivariate distributions, such as the case for matrix-variate Student-t distribution. In this paper, we propose a unified framework which is used not only to introduce a novel multivariate Student-t process regression model (MV-TPR) for multi-output prediction, but also to reformulate the multivariate Gaussian process regression (MV-GPR) that overcomes some limitations of the existing methods. Both MV-GPR and MV-TPR have closed-form expressions for the marginal likelihoods and predictive distributions under this unified framework and thus can adopt the same optimization approaches as used in the conventional GPR. The usefulness of the proposed methods is illustrated through several simulated and real-data examples. In particular, we verify empirically that MV-TPR has superiority for the datasets considered, including air quality prediction and bike rent prediction. At last, the proposed methods are shown to produce profitable investment strategies in the stock markets.


Introduction
Over the last few decades, Gaussian processes regression (GPR) has been proven to be a powerful and effective method for non-linear regression problems due to many favorable properties, such as simple structure of obtaining and expressing uncertainty in predictions, the capability of capturing a wide variety of behaviour by parameters, and a natural Bayesian interpretation [1,2].In 1996, Neal [3] revealed that many Bayesian regression models based on neural network converge to Gaussian processes (GP) in the limit of an infinite number of hidden units [4].GP has been suggested as a replacement for supervised neural networks in non-linear regression [5,6] and classification [5].Furthermore, GP has excellent capability for forecasting time series [7,8].
Despite the popularity of GPR in various modelling tasks, there still exists a conspicuous imperfection, that is, the majority of GPR models are implemented for single response variables or considered independently for multiple responses variables without consideration of their correlation [2,9].In order to resolve the multi-output prediction problem, Gaussian process regression for vector-valued function is proposed and regarded as a pragmatic and straightforward method.The core of this method is to vectorize the multiresponse variables and construct a "big" covariance, which describes the correlations between the inputs as well as between the outputs [2,9,10,11].This modelling strategy is feasible due to the fact that the matrix-variate Gaussian distributions can be re-formulated as multivariate Gaussian distributions [10,12].Intrinsically, Gaussian process regression for vector-valued function is still a conventional Gaussian process regression model since it merely vectorizes multi-response variables, which are assumed to follow a developed case of GP with a reproduced kernel.As an extension, it is natural to consider more general elliptical processes models for multi-output prediction.However, the vectorization method cannot be used to extend multi-output GPR because the equivalence between vectorized matrix-variate and multivariate distributions only exists in Gaussian cases [12].
To overcome this drawback, in this paper we propose a unified framework which: (1) is used to introduce a novel multivariate Student−t process regression model (MV-TPR) for multi-output prediction, (2) is used to reformulate the multivariate Gaussian process regression (MV-GPR) that overcomes some limitations of the existing methods, (3) can be used to derive regression models of general elliptical processes.Both MV-GPR and MV-TPR have closed-form expressions for the marginal likelihoods and predictive distributions under this unified framework and thus can adopt the same optimization approaches as used in the conventional GPR.The usefulness of the proposed methods is illustrated through several simulated examples.Furthermore, we also verify empirically that MV-TPR has superiority in the prediction based on some widely-used datasets, including air quality prediction and bike rent prediction.The proposed methods are then applied to stock market modelling which shows that the profitable stock investment strategies can be obtained.
The rest of the paper is organised as follows.Section 2 introduces some preliminaries of matrix-variate Gaussian and Student−t distributions with their useful properties.Section 3 presents the unified framework to reformulate the multivariate Gaussian process regression and to derive the new multivariate Student−t process regression models.Some numerical experiments by the simulated data and real data and the applications to stock market investment are presented in Section 4. Conclusion and discussion are given in Section 5.

Backgrounds and notations
Matrix-variate Gaussian and Student−t distributions have many useful properties, as discussed in the literatures [12,13,14].For completeness and easy referencing, below we list some of them which will be used in this paper.

Matrix-variate Gaussian distribution
Definition 2. 1 The random matrix X ∈ R n×d is said to have a matrix-variate Gaussian distribution with mean matrix M ∈ R n×d and covariance matrix Σ ∈ R n×n , Ω ∈ R d×d if and only if its probability density function is given by p(X|M, Σ, Ω) = (2π) − dn 2 det(Σ) where etr(•) is exponential of matrix trace, Ω and Σ are positive semi-definite.It is denoted X ∼ MN n,d (M, Σ, Ω).Without loss of clarity, it is denoted X ∼ MN (M, Σ, Ω).
Like multivariate Gaussian distribution, matrix-variate Gaussian distribution also possesses several important properties as follows.
The matrix-variate Gaussian is related to the multivariate Gaussian in the following way.
where vec(•) is vector operator and ⊗ is Kronecker product (or called tensor product).
Furthermore, the matrix-variate Gaussian distribution is consistent under the marginalization and conditional distribution.

Matrix-variate Student−t distribution
where Ω and Σ are positive semi-definite, and ).

Theorem 2.4 (Expectation and covariance)
Remark 2.1 It can be seen that matrix-variate Student−t distribution has many properties similar to matrix-variate Gaussian distribution, and it converges to matrix-variate Gaussian distribution if its degree of freedom tends to infinity.However, matrix-variate Student−t distribution lacks the property of vectorizability (Theorem 2.2) [12].As a consequence Student−t process regression for multiple outputs cannot be derived by vectorizing the multi-response variables.In the next section, we propose a new framework to introduce multivariate Student−t process regression model.

Multivariate Gaussian process regression (MV-GPR)
If f is a multivariate Gaussian process on X with vector-valued mean function u : X → R d , covariance function (also called kernel) k : X × X → R and positive semi-definite parameter matrix Ω ∈ R d×d , then any finite collection of vector-valued variables have a joint matrix-variate Gaussian distribution, where f , u ∈ R d are row vectors whose components are the functions { f i } d i=1 and {µ i } d i=1 respectively.Furthermore, M ∈ R n×d with M ij = µ j (x i ), and Σ ∈ R n×n with Σ ij = k(x i , x j ).Sometimes Σ is called column covariance matrix while Ω is row covariance matrix.We denote f ∼ MGP (u, k, Ω).
In conventional GPR methods, the noisy model y = f (x) + ε is usually considered.However, for Student−t process regression such a model is analytically intractable [16].Therefore we adopt the method used in [16] and consider the noise-free regression model where the noise term is incorporated into the kernel function.
Given n pairs of observations {(x i , y i )} n i=1 , x i ∈ R p , y i ∈ R 1×d , we assume the following model where and Note that the second term in (3) represents the random noises.
We assume u = 0 as commonly done in GPR.By the definition of multivariate Gaussian process, it yields that the collection of functions [f (x 1 ), . . ., f (x n )] follow a matrix-variate Gaussian distribution To predict a new variable T and the predictive targets f * are given by where Thus, taking advantage of conditional distribution of multivariate Gaussian process, the predictive distribution is where Additionally, the expectation and the covariance are obtained,

Kernel
Although there are two covariance matrices in the above regression model, the column covariance and the row covariance, only the column covariance depends on inputs and is considered as kernel since it contains our presumptions about the function we wish to learn and define the closeness and similarity between data points [17].As in conventional GPR, the choice of kernels has a profound impact on the performance of multivariate Gaussian process regression (as well as multivariate Student−t process regression introduced later).A wide range of useful kernels have been proposed in the literature, such as linear, rational quadratic, and Matérn [17].But the Squared Exponential (SE) kernel is the most commonly used due to its simple form and many desirable properties such as smoothness and integrability with other functions, although it could oversmooth the data, especially financial data.
The Squared Exponential (SE) kernel is defined as where s 2 f is the signal variance and can also be considered as an output-scale amplitude and the parameter ℓ is the input (length or time) scale [18].The kernel can also be defined by Automatic Relevance Determination (ARD) where Θ is a diagonal matrix with the element components {ℓ 2 i } p i=1 , which represents the length scales for each corresponding input dimension.
For convenience and the purpose of demonstration, SE kernel is used in all our experiments where there is only one input variable whilst SEard is used in those with multiple input variables.It should be noted that there is no technical difficulty to use other kernels in our models.

Parameter estimation
The hyper-parameters involved in the kernel of MV-GPR need to be estimated from the training data.Although Monte Carlo methods can perform GPR without the need of estimating hyperparameters [6,8,19,20], the common approach is to estimate them by means of maximum marginal likelihood due to the high computational cost of Monte Carlo methods.Ω is an extra parameter compared to the conventional GPR model, hence the unknown parameters include the hyper-parameters in the kernel, the noise variance σ 2 n and the row covariance parameter matrix Ω.
Because Ω is positive semi-definite, it can be denoted as Ω = ΦΦ T , where To guarantee the uniqueness of Φ, the diagonal elements are restricted to be positive and denote In MV-GPR model, the observations follow a matrix-variate Gaussian distribution Y ∼ MN n,d (0, K ′ , Ω) where K ′ is the noisy column covariance matrix with element [K ′ ] ij = k ′ (x i , x j ) so that K ′ = K + σ 2  n I where K is noise-free column covariance matrix with element [K] ij = k(x i , x j ).As we know there are hyperparameters in the kernel k so that we can denote K = K θ .The hyper-parameter set denotes According to the matrix-variate distribution, the negative log marginal likelihood of observations is The derivatives of the negative log marginal likelihood with respect to parameter σ 2 n , θ i , φ ij and ϕ ii are as follows where T , E ij is the d × d elementary matrix having unity in the (i,j)-th element and zeros elsewhere, and J ii is the same as E ij but with the unity being replaced by e ϕ ii .The details are provided in A.2. Hence standard gradient-based numerical optimisation techniques, such as Conjugate Gradient method, can be used to minimise the negative log marginal likelihood function to obtain the estimates of the parameters.Note that since the random noise is incorporated into the kernel function the noise variance is estimated alongside the other hyper-parameters.

Comparison with the existing methods
Compared with the existing multi-output GPR methods [11,2,9], our proposed method possesses several advantages.
Firstly, the existing methods have to vectorize the multi-output matrix in order to utilise the GPR models.It is complicated and not always workable if the numbers of outputs and observations are large.In contrast, our proposed MV-GPR has more straightforward form where the model settings, derivations and computations are all directly performed in matrix form.In particular, we use column covariance (kernel) and row covariance to capture all the correlations together in the multivariate outputs, rather than assuming a separate kernel for each output and constructing a "big" covariance matrix by Kronecker product as done in [2].
Secondly, the existing methods rely on the equivalence between vectorized matrix-variate Gaussian distribution and multivariate Gaussian distribution.However, this equivalence does not exist for other elliptical distributions such as matrix-variate Student−t distribution [12].Therefore, the existing methods for multi-output Gaussian process regression cannot be applied to Student−t process regression.On the other hand, our proposed MV-GPR is based on matrix forms directly and does not require vectorization so it can naturally be extended to MV-TPR as we will do in the next subsection.
Therefore, our proposed MV-GPR provides not only a new derivation of the multi-output Gaussian process regression model, but also a unified framework to derive more general elliptical processes models.

Multivariate Student−t process regression (MV-TPR)
In this subsection we propose a new nonlinear regression model for multivariate response, namely multivariate Student−t process regression model (MV-TPR), using the framework discussed in the previous subsections.MV-TPR is an extension to multi-output GPR, as well as an extension to the univariate Student−t process regression proposed in [16].
By definition, if f is a multivariate Student−t process on X with parameter ν > 2, vector-valued mean function u : X → R d , covariance function (also called kernel) k : X × X → R and positive semi-definite parameter matrix Ω ∈ R d×d , then any finite collection of vector-valued variables have a joint matrix-variate where f , u ∈ R d are row vectors whose components are the functions Therefore MV-TPR model can be formulated along the same line as MV-GPR based on the definition of multivariate Student−t process.We present the model briefly below. Given where ν is the degree of freedom of Student−t process and the remaining parameters have the same meaning of MV-GPR model.Consequently, the predictive distribution is obtained as where According to the expectation and the covariance of matrix-variate Student−t distribution, the predictive mean and covariance are given by In the MV-TPR model, the observations are followed by a matrix-variate Student−t distribution Y ∼ MT n,d (ν, 0, K ′ , Ω).The negative log marginal likelihood is Therefore the parameters of MV-TPR contains all the parameters in MV-GPR and one more parameter: the degree of freedom ν.The derivatives of the negative log marginal likelihood with respect to parameter ν,σ 2  n , θ i , φ ij and ϕ ii are as follows where is the derivative of the function ln Γ n (•) with respect to ν.The details are provided in A.3.

Remark 3.1
It is known that the marginal likelihood function in GPR models is not usually convex with respect to the hyper-parameters, therefore the optimization algorithm may converge to a local optimum whereas the global one provides better result [21].As a result, the optimized hyper-parameters obtained by maximum likelihood estimation and the performance of GPR models may depend on the initial values of the optimization algorithm [6,8,22,20].A common strategy adopted by most GPR practitioners is a heuristic method.That is, the optimization is repeated using several initial values generated randomly from a prior distribution based on their expert opinions and experiences, for example, using 10 initial values randomly selected from a uniform distribution.The final estimates of the hyperparameters are the ones with the largest likelihood values after convergence [6,8,22].Further discussion on how to select suitable initial hyper-parameters can be found in [22,23].In our numerical experiments, the same heuristic method is used for both MV-GPR and MV-TPR.

Experiments
In this section, we demonstrate the usefulness of MV-GPR and our proposed MV-TPR using some numerical examples, including simulated data and real data.

Simulated Example
We first use simulation examples to evaluate the quality of the parameter estimation and the prediction performance of the proposed models.

Evaluation of parameter estimation
We generate random samples from a bivariate Gaussian process y ∼ MGP (0, k ′ , Ω), where k ′ is defined as in (3) with the kernel k SE .The hyper-parameters in k SE are set as [ℓ, s 2 f ] = [0.5, 2.5] and Ω = 1 0.8 0.8 1 .The variance of the random noise in k ′ takes values σ 2 n = 0.1, 0.05 and 0.01.As explained in Section 3.1, in our models the random noises are included in the kernel, therefore no additional random errors can be added when the random samples are generated; otherwise two random error terms will result in identifiability issues in the parameter estimation.The covariate x has 100 equally spaced values in [0, 1].We utilize the heuristic method discussed in the Remark 3.1 to estimate the parameters.The experiment is repeated 20 times and we use the parameter Median Relative Error (pMRE) as a measure of the quality of the estimates [24]: where | • | is the absolute value, θi is the parameter estimates in repetition i and θ is the true parameter.The results are shown in Table 1.The similar experiment is also conducted for MV-TPR, where the samples are generated from a bivariate Student−t process y ∼ MT P (ν, 0, k ′ , Ω) with the same true parameters as above and ν = 3.The results are reported in Table 2.
It can be seen that the length scale ℓ 2 is well estimated in both cases, while the estimates for the parameters in the row covariance matrix (ϕ 11 , ϕ 22 and φ 12 ) are not as good but reasonable.This may be because the conjugate gradient optimization algorithm does not manage to reach the global maxima of the likelihood function.Better results may be achieved using optimal design for parameter estimation as discussed in [24] In general the estimates of the parameters are closer to the true values as the noise level decreases, but the improvement in terms of estimation accuracy is not significant.

Evaluation of prediction accuracy
Now we consider a simulated data from two specific functions.The true model used to generate data is given by (1) , ε (2) ], where the random noise [ε (1) , ε (2) ] ∼ MGP (0, k SE , Ω).We select k SE with parameter [ℓ, s For model training, we use fewer points with one part missing so that the zth training data points with z = {3r + 1} 12  r=1 ∪ {3r + 2} 32 r=22 are selected for both y 1 and y 2 .The prediction is then performed for all 100 covariate values in [-10, 10].The RMSEs between the predicted values and the true ones from f 1 (x) and f 2 (x) are calculated.For comparison, the conventional GPR and TPR models are also applied to the two outputs independently.The above experiment is repeated 1000 times and the ARMSE (Average Root Mean Square Error), defined by , is calculated.Here y ij is the jth observation in the ith repetition and ŷij is the corresponding predicted value.The results are reported in Table 3 and an example of prediction is demonstrated in Figure 1.(h) TPR (y 2 ) Figure 1: Predictions for MV-GP noise data using different models.From panels (a) to (d): predictions for y 1 by MV-GPR, MV-TPR, GPR and TPR.From panels (e) to (h): predictions for y 2 by MV-GPR, MV-TPR, GPR and TPR.The solid blue lines are predictions, the solid red lines are the true functions and the circles are the observations.The dash lines represent the 95% confidence intervals The similar experiment is also conducted for the case where the random noise follows a multivariate Student−t process [ε (1) , ε (2) ] ∼ MT P (3, 0, k SE , Ω).We select k SE with parameter [ℓ, s 2 f ] = [1.001,5] and Ω = 1 0.25 0.25 1 .The resulted ARMSEs are presented in Table 4 and an example of prediction is demonstrated in Figure 2.  Figure 2: Predictions for MV-TP noise data using different models.From panels (a) to (d): predictions for y 1 by MV-GPR, MV-TPR, GPR and TPR.From panels (e) to (h): predictions for y 2 by MV-GPR, MV-TPR, GPR and TPR.The solid blue lines are predictions, the solid red lines are the true functions and the circles are the observations.The dash lines represent the 95% confidence intervals From the tables and figures above, it can be seen that the multivariate process regression models are able to discover a more desired pattern in the gap compared with the conventional GPR and TPR models used independently.It also reveals that taking the correlations between the two outputs into consideration improves the accuracy of prediction compared with the methods of modelling each output independently.In particular, MV-TPR performs better than MV-GPR in the predictions for both types of noisy data whilst the performances by TPR and GPR are similar.
It is not surprising that in general MV-TPR works better than MV-GPR when the outputs have dependencies, because the former has more modelling flexibility with one more parameter which captures the degree of freedom.Theoretically, MV-TP converges to MV-GP if the degree of freedom tends to infinity, and to some extent MV-GPR is a special case of MV-TPR.In the above experiment because the observations, even generated from MV-GP, may contain outliers and the sample size (23 training points) is small, MV-TPR with a large degree of freedom may fit the data better than MV-GPR and hence provides a better prediction.These findings coincide with those in [16] in the context of univariate Student−t process regressions.On the other hand, while a training sample of size 23 may be too small for two-dimensional processes such that MV-TPR performs better than MV-GPR, it may be large enough for one-dimensional processes, leading to similar performances by TPR and GPR.
It is noted that the predictive variance by MV-GPR is much smaller than the independent GPR model.This is likely due to the loss of information in the independent model.

Real data examples
We further test our proposed methods on two real datasets1 .The selected mean function is zero-offset and the selected kernel is SEard.Before the experiments are conducted, all the data have been normalized by ỹi = where µ and σ are the sample mean and standard deviation of the data {y i } n i=1 respectively.

Bike rent prediction
This dataset contains the hourly and daily count of rental bikes between years 2011 and 2012 in Capital Bikeshare System with the corresponding weather and seasonal information [25].There are 16 attributes.We test our proposed methods for multi-output prediction based on daily count dataset.After deleting all the points with missing attributes, we use the first 168 data points in the season Autumn because the data is observed on a daily basis (1 week = 7 days) and the whole dataset is divided into 8 subsets (each subset has 3 weeks' data points).In the experiment, the input comprises 8 attributes, including normalized temperature, normalized feeling temperature, normalized humidity, normalized wind speed, whether the day is holiday or not, day of the week, working day or not and weathersit.The output consists of 2 attributes, including the count of casual users (Casual) and the count of registered users (Registered).The cross-validation method is taken as k−fold, where k = 8.Each subset is considered as a test set and the remaining subsets are considered as a training set.Four models, including MV-GPR, MV-TPR, GPR (to predict each output independently) and TPR (to predict each output independently) are applied to the data and predictions are made based on the divided training and test sets.The process is repeated for 8 times, and for each subset's prediction, the MSE (mean square error) and the MAE (mean absolute error) is calculated.The medians of the 8 MSEs and MAEs are then used to evaluate the performance for each output.Finally, the maximum median of all the outputs (MMO) is used to evaluate the overall performance of the multi-dimensional prediction.The results are shown in Table 5.It can be seen that MV-TPR significantly outperforms all the other models in terms of MSE and MAE, and MV-GPR performs the second best, whilst TPR is slightly better than GPR.

Air quality prediction
The dataset contains 9,358 instances of hourly averaged responses from an array of 5 metal oxide chemical sensors embedded in an Air Quality Chemical Multisensor Device with 15 attributes [26].We delete all the points with missing attributes (887 points remaining).The first 864 points are considered in our experiment because the data is hourly observed (1 day = 24 hours) and the whole data set is divided into 9 subsets (each subset has 4-days' data points, totally 864 data points).In the experiment, the input comprises 9 attributes, including time, true hourly averaged concentration CO in mg/m 3 (COGT), true hourly averaged overall Non Metanic HydroCarbons concentration in microg/m 3 (NMHCGT), true hourly averaged Benzene concentration in microg/m 3 (C6H6GT), true hourly averaged NOx concentration in ppb (NOx), true hourly averaged NO2 concentration in microg/m 3 (NO2), absolute humidity (AH), temperature (T) and relative humidity (RH).The output consists of 5 attributes, including PT08.S1 (tin oxide) hourly averaged sensor response, PT08.S2 (titania) hourly averaged sensor response, PT08.S3 (tungsten oxide) hourly averaged sensor response, PT08.S4 (tungsten oxide) hourly averaged sensor response and PT08.S5 (indium oxide) hourly averaged sensor response.The cross-validation method is taken as k−fold, where k = 9.The remaining modelling procedure is the same as in the bike rent prediction experiment, except that k is 9 so that the process is repeated 9 times.The results are shown in Table 6.It can be observed that MV-TPR consistently outperforms MV-GPR, and in overall terms, MV-GPR does not perform as well as independent GPR.These results can likely be explained as follows.In fact, in our proposed framework, all the outputs are considered as a whole and the same kernel (and hyperparameters) is used for all the outputs, which may not be appropriate if different outputs have very different patterns (which is the case in this example).On the other hand, in the independent modelling different GPR models are used for different outputs and they can have different hyperparameters even though the kernel is the same, which may result in better prediction accuracy.It is noted that in this example MV-TPR still works better than MV-GPR because the former offers more modelling flexibility and is thus more robust to model misspecification.This finding is consistent with that by [16] for the univariate Student−t process.

Application to stock market investment
In the previous subsections, the examples show the usefulness of our proposed methods in terms of more accurate prediction.Furthermore, our proposed methods can be applied to produce trading strategies in the stock market investment.It is known that the accurate prediction of future for an equity market is almost impossible.Admittedly, the more realistic idea is to make a strategy based on the Buy&Sell signal in the different prediction models [27].In this paper, we consider a developed Dollar 100 (dD100) as a criterion of the prediction models.The dD100 criterion is able to reflect the theoretical future value of $100 invested at the beginning, and traded according to the signals constructed by predicted value and the reality.The details of dD100 criterion are described in Section 4.4.2.Furthermore, the equity index is an important measurement of the value of a stock market and is used by many investors making trades and scholars studying stock markets.The index is computed from the weighted average of the selected stocks' prices, so it is able to describe how the whole stock market in the consideration performs in a period and thus many trading strategies of a stock or a portfolio have to take the information of the index into account.As a result, our experimental predictions for specific stocks are based on the indices as well.

Data preparation
We obtain daily price data, containing opening, closing, and adjusted closing for the stocks (the details are shown in Section 4.
Inter-day log return: where ACP i is the adjusted closing price of the ith day (i > 1), CP i is the closing price of the ith day, and OP i is the opening price of the ith day.Therefore, there are totally 503 daily log returns and log inter-day returns for all the stocks and indices from 2013 to 2014.

Prediction model and strategy
The sliding windows method is used for our prediction models, including GPR, TPR, MV-GPR, and MV-TPR, based on the indices, INDU, SPX, and NDX.The training sample is set as 303, which is used to forecast for the next 10 days, and the training set is updated by dropping off the earliest 10 days and adding on the latest 10 days when the window is moved.The sliding-forward process was run 20 times, resulting in a total of 200 prediction days, in groups of 10.The updated training set allows all the models and parameters to adapt the dynamic structure of the equity market [27].Specifically, the inputs consist of the log returns of 3 indices, the targets are multiple stocks' log returns and Standard Exponential with automatic relevance determination (SEard) is used as the kernel for all of these prediction models.It is noteworthy that the predicted log returns of stocks are used to produce a buy or sell signal for trading rather than to discover an exact pattern of the future.The signal BS produced by the predicted log returns of the stocks is defined by i=1 are the predicted log returns of a specific stock, {LR i } 200 i=1 are the true log returns while {ILR i } 200 i=1 are the inter-day log returns.The Buy&Sell strategy relying on the signal BS is described in Table 7.

Decision Condition
No action is taken for the rest of the option It is noted that the stocks in our experiment are counted in Dollar rather than the number of shares, which means in theory we can precisely buy or sell a specific Dollar valued stock.For example, if the stock price is $37 when we only have $20, we can still buy $20 valued stock rather than borrowing $17 and then buy 1 share.Furthermore, it is also necessary to explain why we choose the signal BS.By the definition, we rewrite it as ), where {ACP i } 200 i=0 are the last 201 adjusted closing prices for a stock, {CP} 200 i=1 are the last 200 closing prices, and {AOP} 200 i=1 are the adjusted opening prices.If BS i > 0, the predicted closing price should be higher than the adjusted opening price, which means we can obtain the inter-day profit by buying the shares at the opening price2 as long as the signal based on our predictions is accurate.Meanwhile, the opposite manipulation based on BS strategy means that we can avoid the inter-day loss by selling decisively at the opening price.Furthermore, the reasonable transaction fee 0.025% is considered in the experiment since the strategy might trade frequently 3 .As a result, this is a reasonable strategy since we can definitely obtain a profit by buying the shares and cut the loss by selling the shares in time only if our prediction has no serious problem.It is also an executable strategy because the decision is made based on the next day's reality and our prediction models.
At last, BS signal varies in different prediction models so that we denote these Buy&Sell strategies based on MV-GPR, MV-TPR, GPR, and TPR model as MV-GPR strategy, MV-TPR strategy, GPR strategy, and TPR strategy, respectively.

Chinese companies in NASDAQ
In recent years, the "Chinese concepts stock" has received an extensive attention among international investors owing to the fast development of Chinese economy and an increasing number of Chinese firms have been traded in the international stock markets [28].The "Chinese concepts stock" refers to the stock issued by firms whose asset or earning have essential activities in Mainland China.Undoubtedly, all these "Chinese concept stocks" are heavily influenced by the political and economic environment of China together.For this reason, all these stocks have the potential and unneglectable correlation theoretically, which is probably reflected in the movement of stock prices.The performance of multiple targets prediction, which takes the potential relationship into consideration, should be better.Therefore, the first real data example is based on 3 biggest Chinese companies described in Table 8.We apply MV-GPR, MV-TPR, GPR and TPR strategies and the results are demonstrated in Figure 3. Furthermore, Table B1, Table B2 and Table B3 in the appendix summarize the results by period for each stock respectively.In particular, the Buy&Sell signal examples for each stock are shown in Table B4, Table B5 and Table B6 respectively, along with other relevant details.
Secondly, the 4 models, MV-GPR, MV-TPR, GPR and TPR, are applied in the same way as in Section 4.4.3 and the ranking of stock investment performance is listed in Table 10 (the details are summarized in Table C1).On the whole, for each stock, there is no doubt that using Buy&Sell strategy is much better than using Buy&Hold strategy regardless of the industrial sector.Specifically, MV-GPR makes a satisfactory performance overall in the sectors, Industrials, Consumer Services and Financials while MV-TPR has a higher ranking in Health Care in general.Further analysis is considered using industrial sector portfolios, which consists of these grouped stocks by the same weight investment on each stock.For example, the Oil & Gas portfolio investment is $100 with $50 shares CVX and $50 shares XOM while the Technology portfolio investment is $100 with the same $20 investment on each stock in the industrial sector Technology.The diverse industry portfolio investment performance ranking is listed in Table 11 (the details are described in Table C2).Apparently, the Buy&Sell strategies performed better than the Buy&Hold strategies.MV-GPR suits better in three industries, including Consumer Goods, Consumer Services, and Financials, followed by TPR which performed best in Oil&Gas and Industrials.The optimal investment strategy in Health Care is MV-TPR while in Technology industry, using GPR seems to be the most profitable.

Conclusion and discussion
In this paper we have proposed a unified framework for multi-output regression and prediction.Using this framework, we introduced a novel multivariate Student−t process regression model (MV-TPR) and also reformulated the multivariate Gaussian process regression (MV-GPR) which overcomes some limitations of the existing methods.It could also be used to derive regression models of general elliptical processes.Under this framework, the model settings, derivations and computations for both MV-GPR and MV-TPR are all directly performed in matrix form.MV-GPR is a more straightforward method compared to the existing vectorization method and can be implemented in the same way as the conventional GPR.Similar to the existing Gaussian process regression for vector-valued function, our models are also able to learn the correlations between inputs and outputs, but with more convenient and flexible formulations.The proposed MV-TPR also possesses closed-form expressions for the marginal likelihood and the predictive distributions under this unified framework.Thus the same optimization approaches as used in the conventional GPR can be adopted.The proposed methods are also applied to stock market modelling and are shown to have the ability to make a profitable stock investment.For the three "Chinese concept stocks", the Buy&Sell strategies based on the proposed models have more satisfactory performances, compared with the Buy&Hold strategies for the corresponding stocks and three main indices in the US; in particular, the strategy based on MV-TPR has outstanding returns for NetEase among three stocks.When applied to the industrial sectors in Dow 30, the results indicate that the strategies based on MV-GPR have generally considerable performances in Industrials, Consumer Goods, Consumer Services, and Financials sectors, while those based on MV-TPR can make maximum profit in Health Care sector.
It is noted that we used the squared exponential kernel for demonstration in all our experiments.However, it can be expected that other kernels or more complicated kernels may lead to better results, especially in the financial data examples, as the SE kernel may oversmooth data; see [17] for more details on the choice of kernels.In this paper we assume that different outputs are observed at the same covariate values.In practice, different responses may be observed at different locations.These cases are, however, difficult for the proposed framework since all the outputs have to be considered as a matrix in our models, rather than as a vector with adjustable length.Another issue worth noting is that, as discussed in Section 4.3.2, in our proposed framework all the outputs are considered as a whole and the same kernel (and hyperparameters) is used for all the outputs, which may not be appropriate if different outputs have very different patterns such as in the air quality example.Moreover, it is also important to further study how to improve the quality of the parameter estimates, for example using optimal design methods for parameter estimation as discussed in [24].All of these problems are worth further investigation and exploration and will be our future works.

A.1 Matrix derivatives
According to the chain rule of derivatives of matrix, there exists [30]: Letting U = f (X), the derivatives of the function g(U) with respect to X are where X is an n × m matrix.Additionally, there are another two useful formulas of derivative with respect to X.
where X is an n × n matrix, A is a constant m × n matrix and B is a constant n × m matrix.

A.2 Multivariate Gaussian process regression
For a matrix-variate observations where actually Σ = K + σ 2 n I As we know there are several parameters in the kernel k so that we can denote K = K θ .The parameter set denotes Θ = {θ 1 , θ 2 , . ..}. Besides, we denote the parameter matrix Ω = ΦΦ T since Ω is positive semi-definite, where To guarantee the uniqueness of Φ, the diagonal elements are restricted to be positive and denote where E ij is the d × d elementary matrix having unity in the (i,j)-th element and zeros elsewhere, and J ii is the same as E ij but with the unity being replaced by e ϕ ii .The derivatives of the negative log likelihood with respect to σ 2 n , θ i , φ ij and ϕ ii are as follows.The derivative with respect to θ i is where The fourth equality is due to the symmetry of Σ. Due to ∂Σ/∂σ 2 n = I n , the derivative with respect to σ 2 n is: Letting where the third equation is due to the symmetry of Ω.Similarly, the derivative with respect to ϕ ii is ∂ ln det(Ω)

A.3 Multivariate Student−t process regression
The negative log likelihood of observations Y ∼ MT Therefore, the derivative of negative log marginal likelihood with respect to θ i is 2 f ] = [1.001,5] and Ω = 1 0.25 0.25 1 .The covariate x has 100 equally spaced values in [-10, 10] so that a sample of 100 observations for y 1 and y 2 are obtained.
4.3 and Section 4.4.4) and three main indices in the US, Dow Jones Industrial Average (INDU), S&P500 (SPX), and NASDAQ (NDX) from Yahoo Finance in the period of 2013 -2014.The log returns of the adjusted closing price and inter-day log returns are obtained by defining Log return:

Table 1 :
The pMREs of MV-GP samples with different noise levels estimated by MV-GPR

Table 2 :
The pMREs of MV-TP samples with different noise levels estimated by MV-TPR

Table 3 :
The ARMSE by the different models (multivariate Gaussian noisy data)

Table 4 :
The ARMSE by the different models (multivariate Student−t noisy data)

Table 5 :
Bike rent prediction results based on MSEs and MAEs

Table 8 :
Three biggest "Chinese concept" stocks

Table 10 :
Stock investment performance ranking under different strategies

Table 11 :
Industry portfolio investment performance ranking under different strategies The usefulness of the proposed methods is illustrated through several numerical examples.It is empirically demonstrated that MV-TPR has superiority in prediction in these examples, including the simulated examples, air quality prediction and bike rent prediction.

Table C2 :
The detailed industry portfolio investment results under different strategies