extfieldsofstudy:
Computer Science
Mathematics

text:
Multivariate Arrival Times with Recurrent Neural Networks for Personalized Demand Forecasting

Access to a large variety of data across a massive population has made it possible to predict customer purchase patterns and responses to marketing campaigns. In particular, accurate demand forecasts for popular products with frequent repeat purchases are essential since these products are one of the main drivers of profits. However, buyer purchase patterns are extremely diverse and sparse on a per-product level due to population heterogeneity as well as dependence in purchase patterns across product categories. Traditional methods in survival analysis have proven effective in dealing with censored data by assuming parametric distributions on inter-arrival times. Distributional parameters are then fitted, typically in a regression framework. On the other hand, neural-network based models take a non-parametric approach to learn relations from a larger functional class. However, the lack of distributional assumptions make it difficult to model partially observed data. In this paper, we model directly the inter-arrival times as well as the partially observed information at each time step in a survival-based approach using Recurrent Neural Networks (RNN) to model purchase times jointly over several products. Instead of predicting a point estimate for inter-arrival times, the RNN outputs parameters that define a distributional estimate. The loss function is the negative log-likelihood of these parameters given partially observed data. This approach allows one to leverage both fully observed data as well as partial information. By externalizing the censoring problem through a log-likelihood loss function, we show that substantial improvements over state-of-the-art machine learning methods can be achieved. We present experimental results based on two open datasets as well as a study on a real dataset from a large retailer.

Abstract-Access to a large variety of data across a massive population has made it possible to predict customer purchase patterns and responses to marketing campaigns. In particular, accurate demand forecasts for popular products with frequent repeat purchases are essential since these products are one of the main drivers of profits. However, buyer purchase patterns are extremely diverse and sparse on a per-product level due to population heterogeneity as well as dependence in purchase patterns across product categories. Traditional methods in survival analysis have proven effective in dealing with censored data by assuming parametric distributions on inter-arrival times. Distributional parameters are then fitted, typically in a regression framework. On the other hand, neural-network based models take a non-parametric approach to learn relations from a larger functional class. However, the lack of distributional assumptions make it difficult to model partially observed data. In this paper, we model directly the inter-arrival times as well as the partially observed information at each time step in a survivalbased approach using Recurrent Neural Networks (RNN) to model purchase times jointly over several products. Instead of predicting a point estimate for inter-arrival times, the RNN outputs parameters that define a distributional estimate. The loss function is the negative log-likelihood of these parameters given partially observed data. This approach allows one to leverage both fully observed data as well as partial information. By externalizing the censoring problem through a log-likelihood loss function, we show that substantial improvements over state-ofthe-art machine learning methods can be achieved. We present experimental results based on two open datasets as well as a study on a real dataset from a large retailer.
Index Terms-Survival Analysis, Time series analysis, Neural networks, Consumer products, Multivariate statistics, Maximum likelihood modeling, Bayesian network models, Forecasting, Marketing

I. INTRODUCTION
Accurately predicting each customer's behavior for each product is useful in direct marketing efforts which can lead to significant advantages for a retailer by driving increased sales, margin, and return on investment [13]. Of special interest are replenishable products such as regularly consumed food products (e.g. milk) or regularly replenished personal care products (e.g. soap). These products drive store traffic, basket size, and customer loyalty, which are of strategic importance in a highly competitive retail environment [22].
Traditional approaches to this problem defines customers to be "alive" if purchases were made. These models make as-sumptions on the distributions of purchase counts and lifetimes of individual customers as well as population heterogeneity [2], [6], [18], [29]. They can be framed in a Cox Proportional Hazards model where covariates are regressed as a multiplicative factor in a baseline hazard function [23]. Solving the maximum likelihood problem yields optimal distributional estimates that model these behaviors. In a scenario with sparse purchase data, this can be useful since non-purchases can reveal information about whether a customer is likely to make a purchase in the future.
However, these models impose strict assumptions such as independence and stationarity. Covariates are often modeled linearly, further restricting the space of possible functional relations that are possible. While these assumptions were essential for tractability purposes, they can no longer be taken for granted when we wish to model highly-correlated, highdimensional and heterogeneous processes.
More recently, Recurrent Neural Networks (RNNs) were applied in predicting arrival times [1], [14], [31]. They leverage the capacity of RNNs to model sequential data with complex temporal dependencies as well as non-linear associations [9], [28]. However, these models do not deal explicitly with the uncertainty of random arrival times and are not able to properly exploit censored data.
In this work, we integrate a probabilistic approach to model partially observed data with RNNs in a novel sequence-tosequence approach to predict multiple inter-purchase times for each customer. The RNNs predict a sequence of distribution parameters for a random series of "times to next purchase" (i.e. inter-arrival times). At each time step, we either observe complete information when an arrival has occurred or incomplete information in a period of non-arrival. This induces a conditional distribution on the partially observed "inter-arrival time", which allows us to maximize a likelihood function to obtain optimal RNN parameters. The Multivariate Arrival Times Recurrent Neural Network model will be referred to as MAT-RNN in this paper.
The efficacy of this approach is shown through experiments performed on data from several benchmark datasets and a large retailer. We see that our model out-performs other state-ofthe-art machine learning approaches in predicting whether a customer made purchases in the next time period. The results show that MAT-RNN performs better in the ROC-AUC metric [25] (average ROC-AUC over per-product predictions) in 4 out of the 5 categories of products considered. Additionally, results on the benchmark and synthetic datasets show comparable performance increases when compared to traditional survival model techniques and RNNs trained on the usual squared-loss metric. Implementations and comparisons on open datasets will be published on a publicly-accessible repository.

II. RELATED WORK
Machine learning approaches to demand forecasting have gained popularity in recent years. Random Forest models and other ensemble methods in particular have been widely deployed and have enjoyed success with binary predictions due to their scalability to wide datasets, ease of training and regularization strategies [20], [30]. However, tree-based methods are difficult to extend to modeling partial information when dealing with sparsely observed inter-arrival times. An approach dealing specifically with survival is DeepSurv [15], where the output of a neural net is fed into a risk function and used in a Cox Proportional Hazards (CPH) model. A CPH model is a survival model where the actual hazard rate for each individual is assumed to be a multiplicative factor of a common hazard function [17]. Even so, the common roadblock to both of these models is that they are unable to use and model sequentially-dependent covariates and survival statuses.

A. Survival Analysis and RNNs
Recurrent Neural Nets (RNN) have been proven to model data with complex sequential dependencies [9]. Their success hinged on the development of a Long Short Term Memory (LSTM) structure [12] that incorporates gates to recursively update an internal state in order to make sequential pathdependent predictions. A common method in recent years is to train RNNs to make point estimates for time-to-event by minimizing a distance-based metric [1], [14], [31]. The downside is that unobserved arrival times cannot be explicitly accounted for in these models. However, non-arrivals can reveal a significant amount of information and is not currently exploited by these RNN implementations.
Combining survival analysis with RNNs has been explored recently in order to exploit time series covariates. One approach uses RNNs as a feature extraction step that feeds into a survival model [19]. The RNN takes covariates and the sequential outputs are pooled and incorporated as factors in a Cox Proportional Hazard model. The model is then fitted end-toend so that RNN parameters and hazard parameters are jointly trained. Experiments showed that this approach achieves an acceptable result in failure time prediction. However, modeling survival status as developing with covariates concurrently is a more general approach and possibly more realistic as well.
B. Sequence-to-Sequence Approach to Model Arrival Times A recent paper described a framework for using likelihoodbased losses in RNNs to model sequential data in a sequenceto-sequence manner [8]. A naive approach using this method requires that we model arrivals as a point process, where arrival intensity models only the arrival counts per unit time. This cannot be adapted easily to predict purchases by individual customers. A variation of this is WTTE-RNN (Weibull Time-To-Event RNN), which models the inter-arrival times directly by minimizing a likelihood function that utilizes partial information [21] and achieves satisfactory results in predicting univariate arrival processes. However, the theoretical framework proposed makes it hard to extend to modeling multivariate arrival processes. It also assumes a memoryless arrival process, which the author has recognized is unrealistic. Even so, this approach highlights the fact that there can be advantages to viewing this arrival times prediction problem in the context of survival analysis and in using RNNs to predict distributional parameters in a sequence-to-sequence approach. The proposed theoretical framework (MAT-RNN) addresses these problems and shows that it can be extended easily to model Multivariate Arrival Times more generally. Not only can this be used in a survival context where there is only a single arrival time (i.e. death or failure), but can also be used in a multiple arrivals setting (i.e. purchases).

III. A MULTIVARIATE ARRIVAL TIMES RECURRENT NEURAL NET MODEL
The proposed model specifies a Recurrent Neural Net (RNN) to output distributional parameters which represent predictions for the remaining time to arrival. By iterating through time for each customer, the RNN outputs sequential distributional estimates for the remaining time until the next purchase arrival, giving a personalized demand forecast. In this model, for each customer and product pair, every interarrival time is assumed to be a realization of a distinct random variable. It can be dependent on other product purchases as well as on earlier purchases of the same product. These interarrival times are assumed to be independent conditioned on observing a latent state.
At each time step, the random variable observed is not the inter-arrival time, but the remaining time to next arrival, since we have observed partial information (i.e. the lower bound of the next inter-arrival time). This remaining time to next arrival is often referred to as a conditional excess random variable [Section III-A], whose distribution is used to derive the log-likelihood at each time step [Section III-B]. We assume that our observed process follows a conditional independence structure where these conditional excess random variables are assumed to be independent given the internal state of the RNN [Section III-C, III-D]. The loss function is defined to be the negative log-likelihood and the optimal RNN parameters under such a setup generate distributional parameters which are most likely to explain the observed data. Hence, RNN outputs at the end of training period are our best distributional estimates for the remaining time to next purchase.

A. The Conditional Excess / Remaining Lifetime Random Variable
We will denote the random variable representing the remaining time till next arrival conditioned on the current information as Z t . This random variable is not the true inter-arrival time, but is instead a version that is conditioned on observing partial information. In a survival analysis framework where there is only one inter-arrival time (i.e. the time of death), this is the remaining lifetime.
Consider an arrival process where W n is the time of the n-th arrival and let W 0 = 0, the start of training period. Let N (t) be the number of arrivals by time t [Equation 1]. Also let Y n be the inter-arrival time of the n-th arrival, which is the difference between consecutive arrival times [ Equation 2].
In the context of our problem, at a particular time t, the number of arrivals observed is N (t) and we wish to predict the subsequent (i.e. the {N (t) + 1}-th arrival) and its inter-arrival time Y N (t)+1 . Let tse(t) (time-since-event) be the amount of time that has elapsed since the last arrival or start of training period, whichever is smaller [Equation 3]. This represents the censoring information that is available to the RNN at each time t. Define tte(t) (time-to-event) be the amount of time remaining until the next arrival or the end of testing period (τ ), whichever is smaller [ Equation 4].
Consider a sample process with three arrivals where W 1 = 16, W 2 = 28, W 3 = 32, such that Y 1 = 16, Y 2 = 12, Y 3 = 4. Also, N (t) is a piecewise constant function which is 0 for before the first arrival at t = 16 and jumps by 1 at each arrival time W n . We plot tse(t), tte(t) for t until τ = 40, which is the end of the training period [ Figure 1]. The remaining time until next arrival (Z t ) is a conditional random variable that depends only on Y N (t)+1 , which is the inter-arrival time of the subsequent arrival. We can thereby define this random variable given the observed information [Equation 5], which is commonly known as a conditional excess random variable [7].
In our approach, we do not model the distribution of Z t directly. Instead, the RNN predict parameters for Y N (t)+1 and the losses are computed based on the distribution of Z t , which has a distribution induced by Y N (t)+1 based on partial information (i.e. tse(t)) 1 . This additional structure helps to simplify the incomplete information problem.
For example, consider Z = Y − t | Y > t. This random variable Z is conditioned on the fact that Y has been observed to exceed t and we are interested in the excess value (i.e. Y −t).

B. Log-Likelihood Computation
There are two cases to properly define the log-likelihood function. When the next arrival time is observed, the likelihood evaluation is P (Z t ∈ [tte(t), tte(t) + 1]), since inter-arrival times are only discretely observed. However, where the time to next arrival is not observed (i.e. no more subsequent arrivals are observed by end of training), the likelihood evaluation is instead P (Z t > tte(t)), namely the survival function. Therefore at each time t, the random variable Y N (t)+1 which has distribution parametrized by θ t , induces a distribution on Z t . Thus the log-likelihood at each time t can be written as follows [Equation 7].
In survival analysis terminology, we recall that Z t is the remaining lifetime random variable (i.e. time to death conditioned on current age). For the uncensored case, the likelihood is the death probability and in the right-censored case, it is the survival probability.
We plot distributional estimates at two times to illustrate these cases [ Figure 2]. In the uncensored log-likelihood computation, we assume that f3 is the density function of Z 3 , which is the predictive distribution for the remaining time until next arrival at time step 3. Since the next arrival is observed to have occurred at time 6, the remaining time to next arrival is 3, so we evaluate f3 at the value 3. In the censored case, we consider the predictive distribution for remaining time until next arrival at time step 7. We note that the next arrival was not observed by end of training period at time step 10 hence the right tail of Z 7 (i.e. ≥ 3) was used to compute the loglikelihood.
(a) If next arrival is observed then log-likelihood is density (b) If next arrival is unobserved then log-likelihood is survival.

Fig. 2. Log-Likelihood Visualization for Different Censoring Status
It is essential that the class of distributions for Y n is assumed to have differentiable and numerically stable forms for the density and survival functions to exploit the back-propagation algorithm for efficient gradient computation [12]. An example is the Weibull distribution parametrized by scale (λ) and shape (k), which is used in our approach, whose survival function is made up of only exp() and power() functions [ Equation 8]. The details of the conditional excess random variables based on Weibull distributions can be found in the Appendix [Section VIII].
C. An RNN with distributional outputs A particular RNN structure that had proven to be successful in modeling sequential data is Long Short Term Memory (LSTM) [12], although it can be replaced by other RNN implementations. At each time (t), the outputs of the LSTM, which is parametrized by Θ, are passed through an activation function so that they are valid parameters of a distribution function (θ t ). Then, the log-likelihood is computed for each time step (l t ) as defined earlier [ Equation 7]. The computational graph is shown [ Figure 3], where h t is the internal state of the LSTM and X t are the covariates at each time t.
The loss under such a set-up is the negative of the loglikelihood. Optimal parameters for the LSTM (Θ) are ones that output a series of distributional estimates θ t that best "explain" the sequence of data observed. In the event of an uncensored arrival time at time t, the optimal choice of weights θ t is one that generates a density that has a peak close to the actual arrival time.  hidden state We can therefore factor the joint distribution of {Z t }, giving the log-likelihood computation for the entire time series as a sum of log-likelihoods at each time, such that we obtain the sum below, for arbitrary events E t . Since E t is determined by the censoring status [ Assuming that the LSTM/RNN is parametrized by Θ, we note that there exists a function g that recursively maps X t to (θ t , h t ) that depends only on Θ [Equation 10]. By substituting h t−1 , we can write l t (θ t ) = l t (g t (Θ)) where g t depends only on Θ, g, {X l } l≤t . Then since the overall log-likelihood is a sum of l t (θ t ), it can be written as a function of only the RNN parameters (Θ) and observed data. The structure of the RNN and the back-propagation algorithm allows us to compute gradients of any order efficiently and therefore find theΘ, the minimizer of the overall observed loss [12].

E. Activation Functions for Distributional Parameters
The neural net outputs must be transformed such that they are parameters of a distribution. In our case, we used a Weibull distribution, which is parametrized by shape and scale parameters, both of which are positive values. We initialized the RNN output for scale at the maximum-likelihood estimate (MLE) for the scale parameter of a Weibull distribution whose shape parameter is 1 as this was found to be useful in preventing likelihood-evaluation errors [21]. We also chose a maximum shape parameter (set at 10) and pass the RNN output for shape through a sigmoid function, which is rescaled and shifted such that σ * : R → (0, 10) and σ * (0) = 1. For the scale parameter, an exponential function is used, which is rescaled such that it maps 0 to the average inter-arrival-time.

F. Extension to Multi-Variate Waiting Times
To model multivariate arrivals, we assume there are p different arrival processes of interest. For the i-th waiting time of interest, we define W i,n to be the time of the n-th arrival of this type and N i (t), Y i,n be likewise defined. We also need to define tse(i, t), tte(i, t) to be that for the i-th type as well. Similarly define its associated conditional excess random variable [ Equation 11].
Then we write the log-likelihoods for each event type where l i,t (θ i,t ) = log P (Z i,t = tte(i, t)) or l i,t (θ i,t ) = log P (Z i,t > tte(i, t)), recalling that the former is for the case where the next arrival is observed while the latter is for the case where the no arrivals are observed until the end of training.
The earlier Bayesian Network structure [ Figure 4] requires minimal modifications as we merely require that the emissions are conditionally independent given h t . This then allows us to compute the log-likelihood at each time as a sum, l t (θ t ) = i l i,t (θ i,t ). Since the LSTM network is still parameterized by Θ, the remaining operations are exactly the same as earlier.
In this way, temporal dependence as well as dependence between the p arrival processes can be modeled by the RNN, whose weights Θ will then be optimized by training data. This allows us to model other outputs as well by appending [K 1,t , . . . , K p,t ] to Z t where K j,t is some other variable of interest for process j at time t.

G. Masking of Non-Inter-Arrival-Times
In multi-variate purchase arrival times, we found that masking sequences observed before the first arrival of each product is useful in preventing numerical errors encountered in stochastic gradient descent. As such, log-likelihoods computed for time steps before the earliest arrival are masked. This ensures that RNN parameters are not updated due to losses incurred during these times.

H. Predicting Time to Next Arrival
Predictions are simple to compute since at each time t, we can use the estimated parameter θ t to compute the expectation of any function of Z t , assuming that Z t is distributed according to θ t . Since we are concerned with the next arrival time after the end of training period (time τ ), we can compute many different values of interest.
For example, if we want to find the predicted probability that the next arrival will occur within γ time after end of training, we can compute P (Z τ ≤ γ). We can also compute a deferred arrival probability, which is the probability that the next arrival will occur within an interval between γ 1 and γ 1 +γ 2 time after end of training given that we know it will not occur within γ 1 time after the end of training. This can be found by computing P (Z τ ∈ [γ 1 , γ 1 + γ 2 ] | Z τ > γ 1 ). The quantities of interest may not necessarily be limited to probabilities (e.g. mean, quantiles of the predictive distribution) and can be extended to generate other analytics for revenue analysis or forecasting that depends on the subsequent purchase time.

IV. EXPERIMENTAL RESULTS
We ran experiments to show the efficacy of the proposed model. Since our model provides distributional estimates for time to next event, we can evaluate its performance in a binary classification framework as well as in a regression context. This is explored by running experiments on two open datasets. The first is store purchases so there are multiple arrivals but with no covariates. The second is failure time prediction with many covariates. Finally, we apply the model to a real dataset of retail purchases in order to compare realworld performances of the considered models. We were unable to find open data for multivariate arrival times with a large number of arrival types but the publicly available code can be easily extended to run on such data.

A. Performance Evaluation
Binary classification performance is evaluated using the Receiver Operating Characteristics Area Under Curve (ROC-AUC), which is equivalent to the Concordance Index (C-Index) commonly used in survival analysis [11]. Point prediction performance is evaluated using specified distance metrics.

B. Recurrent Neural Net Structure
For all experiments in this section, a simple structure for our neural net is chosen where there are only three stacked layers, with two LSTM layers of size W followed by a densely connected layer of size 2p, where p is the number of arrival processes. The densely connected layer transforms the LSTM outputs to a vector of length 2p. We chose a Weibull distribution as the family of distributions for our inter-arrival times, which has two parameters, namely the scale and shape. In MAT-RNN, the outputs are then passed through an activation layer [Section III-E]. For squared-loss RNNs, the activation is passed through a softplus layer since time to arrivals are nonnegative. A masking layer is applied prior to the other layers on time indices prior to the initialization of the time series so that losses incurred during those time steps do not affect optimization.

Fig. 5. Recurrent Network Setup for a p-variate Arrival Process with W -long Hidden State
This structure is the same for other neural network based models used for benchmark comparison. The networks are trained with Adam [16] and non-default training parameters are specified where possible. Gradients are component-wise clipped at 5 to prevent numerical issues. Implementation was done using Keras [4] as our TensorFlow [5] wrapper. Evaluation metrics and ensemble predictors are implemented using scikit-learn [24]. Implementations and comparisons on open datasets will be published on a publiclyaccessible repository 2 .

C. Comparisons on Open Datasets
We show the flexibility of the proposed MAT-RNN model with two open dataset benchmarks. These two problems are often tackled with different models since the prediction problem is different. The proposed MAT-RNN model however, can be adapted to solve these problems since they can be modeled by a distributional approach to inter-arrival times.
The first example is the CDNOW dataset [Section IV-C1]. We consider a binary classification problem where we predict if purchases are made during a testing period. Predictions from MAT-RNN can be computed as the probability that the inter-arrival time occurs before the end of the testing period. Data available is limited to only the transaction history (i.e. purchase date, purchase quantity, customer id, etc) without other covariate data. This type of data is particularly suited to the simple Pareto/NBD type models discussed earlier. We show that MAT-RNN out-performs on a dataset even with no covariates.
The second example is based on the CMAPSS dataset [Section IV-C2] where we predict the remaining useful lifetime, or the time to failure. Predictions from MAT-RNN can be computed as the mode, mean or some other function of the inter-arrival time distribution. The training data is an uncensored time series where sensor readings and operational settings are collected until the engine fails. A customized loss function was used to evaluate models in the PHM08 competition [26], which we will include for our evaluation. Since the training data is fully observed, we would expect the RNN model to perform well.
1) Binary Classification Comparison on CDNOW: Purchase transactions are available from the CDNOW dataset 3 [10] where number of customer purchases are recorded. Only transaction dates, purchase counts and transaction value are available as covariates. Our model is trained on weekly purchases and hidden layers are 1-long (i.e. W = 1 [Section IV-B]). Binary classification performance is compared to the Pareto/NBD model, which is a classical demand forecasting model [29] using the lifetimes package [3]. This dataset is often used as an example where Pareto/NBD type models do well since there's limited covariate data available and there's only a single event type.
With W = 1, there are only 32 trainable parameters in the MAT-RNN model. The training period is set at 1.5 years, from 1997-01-01 to 1998-05-31. Predictions are made for customer purchases within a month of the end of training (i.e. before 1998-06-30). We see that MAT-RNN achieves an ROC-AUC of 0.84 which is better compared to 0.80 that is obtained using the Pareto/NBD estimate for the "alive" probability [ Figure 6]. Results from WTTE-RNN are similar and yields almost the same ROC-AUC curve. This is likely due to the mostly memoryless nature of these customer purchases. It can be seen that this approach of integrating a survivalbased maximum log-likelihood method with a RNN can yield improved prediction accuracy even with a small number of weights and on a small dataset.  [27], is a high dimensional dataset on engine performance with 26 sensor measurements and operational settings. In the training model, the engines are run until failure. In the testing model, data is recorded until a time prior to failure. The goal is to predict the remaining useful life (RUL) for these engines. We used the first set of engine simulations in the dataset (train_FD001.txt) which has 100 uncensored time series of engines that were run until failure. The maximum cycles run before failure was found to be at 363. Time series for each engine was segmented into sliding windows of window length 78, resulting in 323 windowed time series each of length 78. For the testing dataset, RNN is run on time series 78 cycles before end of observation.
A custom loss function was defined for the PHM08 conference competition that was based on this dataset, where overestimation is more heavily penalized [26]. The loss function is defined as follows, where d is the predicted RUL subtracted by the true RUL [ Equation 12].
We compare the performances of the RNN models since they are able to make point predictions fro RUL. The SQ-LOSS model was used with a softplus activation layer scaled by average failure time and the weights are fitted on squared loss. The RNN models were trained with W = 64 [Section IV-B], where there are 2 hidden layers each of Wlong. A variety of training iterations and learning rates were used, and a grid search was performed over all combinations. Performance is evaluated on the test set based on the root mean squared loss metric (rMSE) as well as the mean custom loss metric (MCL) [Equation 12]. Testing set losses are presented in a table [Table IV-C2]. We find that MAT-RNN is generally easier to train, achieving the smallest losses under almost all hyper-parameters other than a slow learning rate and small training iteration. We visualize the predictions of the best performing model, achieved by MAT-RNN with iters = 1e4, lr = 1e − 5 [ Figure  7]. As expected, we find that the predicted density of the remaining lifetime to be highest at the true RUL. Taking point estimate as the model of the predicted density, we can then plot predicted RUL against true RUL, which indicates that we have a good predictive accuracy.

D. Demand Forecasts for a Large Retailer
We investigate the predictive performance for purchases of of a few baskets of goods sold by a large retailer. The time resolution of our dataset is on a weekly level. Training data is available over roughly 1.5 years, which gives us 78 weeks of training data from 2014-01-01 to 2015-06-30. Performance of our MAT-RNN model is measured on the binary classification problem of predicting whether a customer purchases the product within 4 weeks after the end of training period from 2015-06-30 to 2015-07-31.
1) Benchmark Models: Even though our model can predict many different quantities of interest [Section III-H], we chose to compare the predictive performance to a few benchmark models in terms of whether an event will arrive within γ time after the end of the training period. These models are namely the Squared-Loss RNN (SQ-RNN) and a Random Forest Predictor (RNG-F), which are state-of-the-art models in personalized demand forecasting [20] and forecasting event arrival times [1]. Models are trained on all customers who bought an item in the basket during the training period and performance is evaluated on this group of customers during the testing period.
• RNG-F is trained by splitting the training period into two periods. Covariates at the end of the first period are fed into the model, which is trained to predict whether subjects purchase in the second period of length γ. A different RNG-F model is trained for each product, but is fed covariate datasets for all products. • SQ-LOSS is trained by setting the loss function as the squared difference between the predicted time-to-arrival and the actual time-to-arrival. An activation function of softplus is applied. Predictions of SQ-LOSS are then compared to the testing period length of γ. If by the end of the training period, SQ-LOSS predicts the next timeto-arrival as s, then the prediction metric is γ−s. For time periods where no actual time-to-arrival was observed (i.e. no further purchases were observed by end of training), loss is set to 0. 2) Covariates: For each customer, at each time period, we compute the Recency, Frequency and Monetary (RFM) metrics which are commonly used in demand modeling [20] at three different levels: namely for all products, in-basket products and each individual product. Recency is the time since last purchase, frequency is the number of repeat purchases and monetary is the amount of money spent on all purchases to date. Included in the covariates are time-since-event (tse(t)) and indicators for whether a first purchase has occurred (pch(t)). We also compute the time-to-event (tse(t)) as well as the censoring status of the next arrival (unc(t)), which are only passed to the loss function. On a per-product level, the types of covariates are limited to only RFM metrics (3 covariates) and transformations of purchase history (2 series). RFM metrics on the category and overall purchase history levels are available as well, but these account for an additional 6 covariates that are shared across the various purchase arrival processes. The total number of covariates for each product is thus 11, of which 6 are shared with other products.
3) Data Summaries for Product Baskets: We selected 5 baskets of popular replenishable products. These are selected from products ranked by a score, where N unique is the number of unique customers and X is the average purchases per customer. score = X * log N unique (13) The selected baskets are bars, deli, floss, pads, soda. Their data summaries are presented [Table II], where µ overall is the average in-basket purchase counts, µ per-sku is the mean over the per-product average purchase counts and p others is the mean over the per-product proportion of buyers who bought another product in-basket. Also note that p trial is the mean over the per-product proportion of trial customers (i.e. those who have made only a single purchase). We note that pads has the highest proportion of trial customers along with the smallest proportion of customers who bought another item in the basket. On the other hand, we find that µ per-sku is roughly median in the baskets considered. This is similar for floss as well. For these categories, it would be reasonable to expect product purchases are strongly dependent. A good joint-prediction model should separate trial purchasers from repeat purchasers who decided to stick to one product after trying another.

4) Performance of Joint Predictions:
Performance is measured based on the ROC-AUC metric where each of the models (i.e. RNG-F, SQ-LOSS, MAT-RNN) predict whether customers who made in-basket purchases will make another in-basket purchase in a 4 week period after the end of a training period of 78 weeks. Framing this as a binary prediction problem allows us to compare these state-of-theart models. The RNN-based models share the same network specifications with W = 36 [Section IV-B] and predict arrival times jointly over different products for each customer. These models are run for 100 iterations with a learning rate set to 1e − 3. The RNG-F model is trained with 100 trees with covariates at week 74 and purchases between week 74 and 78 but predicts purchases for only one product at a time using default parameters in scikit-learn. As such, a separate RNG-F model is trained for each product. We note that RNG-F model is favored since the training metric is directly related to the testing metric, whereas the RNN-based models train and predict based on arrival times.
We are interested in the performance of each model for every product in the basket so there are multiple ROC-AUC metrics. The results are presented in terms of summary statistics for ROC-AUCs for each item in the basket [ Table  III]. In our testing, we found that the MAT-RNN model almost always dominates in the ROC-AUC metric for every category other than bars and deli, which has the smallest number of customers. Even so, MAT-RNN still performs the best in terms of average ROC-AUC among products in each category [ Table III] other than bars.
The number of products for which ROC-AUC has improved over RNG-F is substantial for MAT-RNN. Excluding bars where only 2 out of 6 products saw improved performance, other categories saw ROC-AUC improvements in more than 60% of the products in-category, with soda and pads showing improvements in all products. The ability to model sequential data and sequential dependence separates MAT-RNN model from RNG-F. Even though RNG-F is trained on the evaluation metric, we find that MAT-RNN almost always performs better in this binary classification task.
Also notable is that the performance difference of MAT-RNN over SQ-LOSS and RNG-F is greatest for the pads category. This is likely due to the large amount of missing data since customers are least likely to buy other products [ Table II]. We also find that SQ-LOSS performs poorly compared to MAT-RNN [Table III], even though these models have the same recurrent structure and are fed the same sequential data. One possible explanation is that the lack of ground truth data has a significant impact on the ability of SQ-LOSS to learn. In cases where event arrivals are sparse or where inter-purchase periods are long, the censored nature of the data gives no ground truth to train SQ-LOSS on. Therefore, even though the recurrent structure makes it possible to model sequential dependence, the structure that MAT-RNN imposes on the problem makes it much easier to make predictions with censored observations. We should note that RNG-F out-performs MAT-RNN for 4 out of 12 deli products and 4 out of 6 bars products. A possible reason is that these categories have a smaller sample size. One limitation of the RNN models is that regularization options are more limited compared to ensemble methods used in RNG-F. Based on the results, it appears that MAT-RNN performs better for more customers [ Table III].

5) Joint Predictions and Individual Predictions:
It is clear that joint predictions enjoy some advantages over individual predictions. As stated earlier, we expect correlations in product purchases to be modeled better through joint modeling. If network structure is the same, then the amount of time required to train a separate model for each product scales linearly with the number of products. The number of parameters in a collection of individual models is also significantly larger than that of a joint model.
We study the advantages of training a joint MAT-RNN model over a collection of individual ones by comparing ROC-AUC performance in the soda basket. The per-product individual models are given the same covariates but trained only on the purchase arrivals of that particular product. The network structure is the same with W = 36, but the final densely connected layer outputs only a vector of size 2, since distributional parameters for one product is required. However, since the collection of single models have different weights for their RNNs, they have approximately 8 times the number of

V. CONCLUSION
We described a method to incorporate a survival analysis approach with recurrent neural nets (RNN) to forecast joint arrival times till next purchase for each individual customer over multiple products. This was achieved by transforming the arrival time problem into a likelihood-maximization one with loose distributional assumptions regarding inter-arrival times. Our experiments show that this leads to significant improvement over current state-of-the-art methods is possible. This is the result of being able to model purchases jointly as well as combining a parametric approach to modeling partially observed information with recent advances in neural networks.

VI. FUTURE WORK
Future work includes modeling more features, as the structure allows for the observed processes to be independent given the parameter generation process. An extra feature to model is first-arrival times, which can help predict whether a customer might make their first purchase in response to marketing campaigns. However, to model this effectively, we require a model-based solution to address the sparsity of observations of first-arrival times.

A. Scaling to Inventory-Wide Predictions
Including thousands of products is infeasible in the current implementation. It would be useful to extend the model to make joint inventory-wide predictions for purchase arrival times. Inventory-wide predictions also brings about other technical problems such as memory usage and ETL (data extraction, transformation, loading). The computation of timesince-event and time-to-event as well as other features to be used in log-likelihood computation [Section III-B] is also problematic in that it transforms a sparse dataset (i.e. transactional information) into a dense dataset (i.e. temporal data). An inventory-wide extension needs to address these issues by either implementing an efficient method to compute these matrices or avoid the problem through a variant of MAT-RNN. Nonetheless, we believe that the current implementation can suffice for the purposes of targeted advertising where we are only interested in predicting demand and pushing sales for a few product lines.

VII. ACKNOWLEDGMENTS
We would like to thank Kanchana Padmanabhan and Erin Wilder at Rubikloud Technologies Inc for their expertise in editing this paper as well as the Data Science and Engineering team for their technical support. This would also not have been possible without the supervision and advice of Professor Nancy Reid and Professor Andrei Badescu at the University of Toronto. Lastly, we'd like to thank the NSERC (Natural Sciences and Engineering Research Council) Engage and Mitacs Accelerate programs for funding this project.
The conditional excess random variable, given that it exceeds s, is W = Y − s|Y > s. Recall the definition of conditional probability in terms of some continuous random variable X 1 , X 2 , for any measurable set A 1 , A 2 , given P (X 2 ∈ A 2 ) > 0: We can therefore derive the conditional excess survival function: S W (t) = P (W > t) = P (Y > s + t|Y > s) = S(s + t)/S(s) = exp −((s + t)/λ) k + (s/λ) k .


text:
Vector and Matrix Optimal Mass Transport: Theory, Algorithm, and Applications

In many applications such as color image processing, data has more than one piece of information associated with each spatial coordinate, and in such cases the classical optimal mass transport (OMT) must be generalized to handle vector-valued or matrix-valued densities. In this paper, we discuss the vector and matrix optimal mass transport and present three contributions. We first present a rigorous mathematical formulation for these setups and provide analytical results including existence of solutions and strong duality. Next, we present a simple, scalable, and parallelizable methods to solve the vector and matrix-OMT problems. Finally, we implement the proposed methods on a CUDA GPU and present experiments and applications.

However, in many applications such as color image processing, there is more than one piece of information associated with each spatial coordinate, and such data can be interpreted as vector-valued or matrix-valued densities. As the classical optimal mass transport works with scalar probability densities, such applications require a new notion of mass transport. For this purpose, Chen et al. [18,19,14,12] recently developed a framework for vector-valued and matrix-valued optimal mass transport. See also [41,40,42,25,45,55] for other different frameworks. For vector-valued OMT, potential applications include color image processing, multi-modality medical imaging, and image processing involving textures. For matrix-valued OMT, we have diffusion tensor imaging, multivariate spectral analysis, and stress tensor analysis.
Several mathematical aspects of vector and matrix-valued OMT were not addressed in the previous work [18,19,14]. As the first contribution of this paper, we present duality and existence results of the continuous vector and matrix-valued OMT problems along with rigorous problem formulations.
Although the classical theory of OMT is very rich, only recently has there been much attention to numerical methods to compute the OMT. Several recent work proposed algorithms to solve the L 2 OMT [3,22,7,28,6,13,17,27] and the L 1 OMT [34,50]. As the second contribution of this paper, we present first-order primal-dual methods to solve the vector and matrix-valued OMT problems. The methods simultaneously solve for both the primal and dual solutions (hence a primal-dual method) and are scalable as they are first-order methods. We also discuss the convergence of the methods.
As the third contribution of this paper, we implement the proposed method on a CUDA GPU and present several applications. The proposed algorithms' simple structure allows us to effectively utilize the computing capability of the CUDA architecture, and we demonstrate this through our experiments. We furthermore release the code for scientific reproducibility.
The rest of the paper is structured as follows. In Section 2 we give a quick review of the classic OMT theory, which allows us to present the later sections in an analogous manner and thereby outline the similarities and differences. In Section 3 and Section 4 we present the vector and matrix-valued OMT problems and state a few theoretical results. In Section 5, we present and prove the analytical results. In Section 7, we present the algorithm. In Section 8, we present the experiments and applications.

Scalar optimal mass transport
The theory of optimal transport [26,53,54] remarkably points out that (2) is equivalent to the following flux minimization problem: where u = (u 1 , . . . , u d ) : Ω → R d is the optimization variable and div x denote the (spatial) divergence operator. Although (3) and (2) are mathematically equivalent, (3) is much more computationally effective as its optimization variable u is much smaller when discretized. It is worth mentioning that OMT in formulation (3) is very close to the problems in compressed sensing. Its objective function is homogeneous degree one and the constraint is linear. It can be observed that for characterizing the OMT, the gradient operator in (3) and divergence operator in (3) play the key roles. Later on, we extend the definition of L 1 OMT problem by extending these differential operators into a general meaning.
The optimization problem (3) has the following dual problem: where the optimization and φ : Ω → R is a function. We write · * for the dual norm of · . It is well-known that, strong duality holds between (3) and (4) in the sense that the minimized and maximized objective values are equal [53]. Therefore, we take either (3) or (4) as the definition of S.
Rigorous definitions of the optimization problems (3) or (4) are somewhat technical. We skip this discussion, as scalar optimal mass transport is standard. In Section 5, we rigorously discuss the vector-OMT problems, so any rigorous discussion of the scalar-OMT problems can be inferred as a special case. To make c j the cost of traversing edge j, the edge weight is defined to be 1/c 2 j for j = 1, . . . , .

Theoretical properties
The algorithm we present in Section 7 is a primal-dual algorithm and, as such, finds solutions to both Problems (3) and (4). This is well-defined as both the primal and dual problems have solutions [53]. Write R + for the set of nonnegative real numbers. Write P(Ω, R) for the space of nonnegative densities supported on Ω with unit mass. We can use S(λ 0 , λ 1 ) as a distance measure between λ 0 , λ 1 ∈ P(Ω, R). The value S : P(Ω, R) × P(Ω, R) → R + defines a metric on P(Ω, R) [53].

Vector optimal mass transport
Next we discuss the vector-valued optimal transport, proposed recently in [19]. The basic idea is to combine scalar optimal mass transport with network flow problems [1].

Gradient and divergence on graphs
Consider a connected, positively weighted, undirected graph G with k nodes and edges. To define an incidence matrix for G, we say an edge {i, j} points from i to j, i.e., i → j, if i < j. This choice is arbitrary and does not affect the final result. With this edge orientation, the incidence matrix D ∈ R k× is For example, the incidence matrix of the graph of Figure 1 is Write ∆ G = −D diag{1/c 2 1 , · · · , 1/c 2 }D T for the (negative) graph Laplacian, where 1/c 2 1 , · · · , 1/c 2 are the edge weights. The edge weights are defined so that and c j represents the cost of traversing edge j for j = 1, . . . , .
We define the gradient operator on G as ∇ G x = diag{1/c 1 , · · · , 1/c }D T x and the divergence operator as div G y = −D diag{1/c 1 , · · · , 1/c }y. So the Laplacian can be rewritten as ∆ G = div G ∇ G . Note that div G = −∇ * G , where ∇ * G is the adjoint of ∇ G . This is in analogy with the usual spatial gradient and divergence operators [36,20,21,19].

Vector optimal mass transport
We say λ : Ω → R k + is a nonnegative vector-valued density with unit mass if Assume λ 0 and λ 1 are nonnegative vector-valued densities supported on Ω with unit mass. We define the optimal mass transport between vector-valued densities as where u : Ω → R k×d and w : Ω → R are the optimization variables, α > 0 is a parameter, and · u is a norm on R k×d and · w is a norm on R . The parameter α represents the relative importance of the two flux terms u and w. We write We call div x the spatial divergence operator. The zero-flux boundary condition is for i = 1, . . . , k. Note that w has no boundary conditions. The optimization problem (5) has the following dual problem: where the optimization variable φ : Ω → R k is a function. We write · u * and · w * for the dual norms of · u and · w , respectively. As stated in Theorem 1, strong duality holds between (5) and (6) in the sense that the minimized and maximized objective values are equal. Therefore, we take either (5) or (6) as the definition of V . In Section 5, we rigorously define the primal and dual problems and prove Theorem 1.

Theoretical properties
The algorithm we present in Section 7 is a primal-dual algorithm and, as such, finds solutions to both Problems (5) and (6). This is well-defined as both the primal and dual problems have solutions. Write P(Ω, R k ) for the space of nonnegative vector-valued densities supported on Ω with unit mass. We can use V ( λ 0 , λ 1 ) as a distance measure between λ 0 , λ 1 ∈ P(Ω, R k ). Theorem 2. V : P(Ω, R k ) × P(Ω, R k ) → R + defines a metric on P(Ω, R k ).

Quantum gradient operator and matrix optimal mass transport
We closely follow the treatment in [18]. In particular, we define a notion of gradient on the space of Hermitian matrices and its dual, i.e., the (negative) divergence.
Some applications of matrix-OMT, such as diffusion tensor imaging, have real-valued data while some applications, such as multivariate spectral analysis, have complex-valued data [51]. To accommodate the wide range of applications, we develop the matrix-OMT with complex-valued matrices.
Write C, H, and S for the set of k × k complex, Hermitian, and skew-Hermitian matrices respectively. We write H + for the set of k × k positive semidefinite Hermitian matrices, i.e., M ∈ H + if v * M v ≥ 0 for all v ∈ C k . Write tr for the trace, i.e. for any M ∈ H, we have tr(M ) = k i=1 M ii . Write C N for the block-column concatenation of N matrices in C, i.e., Z ∈ C N if Z = [Z * 1 · · · Z * N ] * and Z 1 , . . . , Z N ∈ C. Define H N and S N likewise. For X, Y ∈ C, we use the Hilbert-Schmidt inner product (This is the standard inner product when we view C as the real vector space R 2n 2 .) For X ∈ C, we use the norm

Quantum gradient and divergence operators
We define the gradient operator, given a L = [L 1 , · · · , L ] * ∈ H , as Define the divergence operator as This is in analogy with the usual spatial gradient and divergence operators. Write ∆ L = div L ∇ L .
This notion of gradient and divergence operators is motivated by the Lindblad equation in Quantum mechanics [18]. The choice of L affects ∇ L . There is no standard way of choosing L. A standing assumption throughout, is that the null space of ∇ L , denoted by ker(∇ L ), contains only scalar multiples of the identity matrix I.
Assume Λ 0 and Λ 1 are nonnegative matrix-valued densities supported on Ω with unit mass.
We define the optimal mass transport between matrix-valued densities as where U : Ω → H d and W : Ω → S are the optimization variables, α > 0 is a parameter, and · u is a norm on H d and · w is a norm on S . The parameter α represents the relative importance of the two flux terms U and W. We write We define the spatial divergence as The zero-flux boundary condition is for i, j = 1, . . . , k. Note that W has no boundary conditions [14]. The optimization problem (7) has the following dual problem: where the optimization variable Φ : Ω → H is a function. We write · u * and · w * for the dual norms of · u and · w , respectively. As stated in Theorem 3, strong duality holds between (7) and (8) in the sense that the minimized and maximized objective values are equal. Therefore we take either (7) or (8) as the definition of M .
To avoid repeating the same argument with different notation, we simply point out that the rigorous definitions of the matrix optimal mass transport problems are analogous to those of the vector setup. So the precise definitions of (7) or (8) can be inferred from the discussion of Section 5.

Theoretical properties
The algorithm we present in Section 7 is a primal-dual algorithm and, as such, finds solutions to both Problems (7) and (8). This is well-defined as both the primal and dual problems have solutions.
Theorem 3. The (infinite dimensional) primal and dual optimization problems (7) and (8) have solutions, and their optimal values are the same, i.e., strong duality holds.

Duality proof
In this section, we establish the theoretical results. For notational simplicity, we only prove the results for the vector-OMT primal and dual problems (5) and (6). Analogous results for the matrix-OMT primal and dual problems (7) and (8) follow from the same logic.
Although the classical scalar-OMT literature is very rich, standard techniques for proving scalar-OMT duality do not simply apply to our setup. For example, Villani's proof of strong duality, presented as Theorem 1.3 of [53], relies on and works with the linear optimization formulation (2). However, our vector and matrix-OMT formulations directly generalize the flux formulation (3) and do not have formulations analogous to (2). We need a direct approach to analyze duality between the flux formulation and the dual (with one function variable), and we provide this in this section.
We further assume Ω ⊂ R d has a piecewise smooth boundary. Write Ω • and ∂Ω for the interior and boundary of Ω. For simplicity, assume Ω as full affine dimensions, i.e., Ω • = Ω.
The rigorous form the dual problem (6) is where W 1,∞ (Ω, R k ) is the standard Sobolev space of functions from Ω to R k with bounded weak gradients. That (9) has a solution directly follows from the Arzelà-Ascoli Theorem.
To rigorously define the primal problem (5) requires more definitions, and we do so later as (11).

Fenchel-Rockafellar duality
Let L : X → Y be a continuous linear map between locally convex topological vector spaces X and Y , and let f : In this framework of Fenchel-Rockafellar duality, d ≤ p , i.e., weak duality, holds unconditionally, and this is not difficult to prove. To establish d = p , i.e., strong duality, requires additional assumptions and is more difficult to prove. The following theorem does this with a condition we can use.

Theorem 5.
[Theorem 17 and 18 [48]] If there is an x ∈ X such that f (x) < ∞ and g is bounded above in a neighborhood of Lx,

Spaces
Throughout this section, · 1 , · 2 , . . . denote unspecified finite dimensional norms. As all finite dimensional norms are equivalent, we do not bother to precisely specify which norms they are. Define Then C 1 (Ω, R k ) is a Banach space equipped with the norm · ∞,∞ . Write M(Ω, R k ) for the space of R k -valued signed finite Borel measures on Ω, and define M(Ω, R k×d ) likewise. Write (C(Ω, R k×d )) * , (C 1 (Ω, R k )) * for the topological dual of C(Ω, R k×d ), C 1 (Ω, R k ), respectively. The standard Riesz-Markov theorem tells us that (C(Ω, R k×d )) * = M(Ω, R k×d ).
Fully characterizing (C 1 (Ω, R k )) * is hard, but we do not need to do so. Instead, we only use the following simple fact. Any g ∈ M(Ω, R k ) defines the bounded linear map φ → Ω φ(x), g(dx) for any φ ∈ C 1 (Ω, R k ). In other words, M(Ω, R k×d ) ⊂ (C 1 (Ω, R k )) * with the appropriate identification.

Operators
We redefine ∇ x : C 1 (Ω, R k ) → C(Ω, R k×d ) so that ∇ x φ is the continuous exntension of the usual ∇ x φ to all of Ω. This makes ∇ x a bounded linear operator. Define the dual (adjoint) operator ∇ * for any φ ∈ C 1 (Ω, R k ) and u ∈ M(Ω, R k×d ). Define the ∇ G (which is simply a multiplication by a R k× matrix) as ∇ G : C 1 (Ω, R k ) → C(Ω, R k ). Since C 1 (Ω, R k ) ⊂ C(Ω, R k ), there is nothing wrong with defining the range of ∇ G to be C(Ω, R k ), and this still makes ∇ G a bounded linear operator. Define the dual (adjoint) operator ∇ * G : M(Ω, R k ) → (C 1 (Ω, R k )) * by identifying ∇ * G with the transpose of the matrix that defines ∇ G . Since ∇ * G is simply multiplication by a matrix, we can further say

Zero-flux boundary condition
Let m : Ω → R k×d a smooth function. Then integration by parts tells us that holds for all smooth ϕ : Ω → R k if and only if m(x) satisfies the zero-flux boundary condition, i.e., m(x)n(x) = 0 for all x ∈ ∂Ω, where n(x) denotes the normal vector at x. Here div x denotes the usual (spatial) divergence. To put it differently, ∇ * x = −div x holds when the zero-flux boundary condition holds. We generalize this notion to measures. We say u ∈ M(Ω, R k×d ) satisfies the zero-flux boundary condition in the weak sense if there is a g ∈ M(Ω, R k ) ⊂ (C 1 (Ω, R k )) * such that holds for all φ ∈ C 1 (Ω, R k ). In other words, u satisfies the zero-flux boundary condition if ∇ * x u ∈ M(Ω, R k ) ⊂ (C 1 (Ω, R k )) * . In this case, we write div x ( u) = g and div x ( u) = −∇ * x ( u). This definition is often used in elasticity theory.

Duality
To establish duality, we view the dual problem (6) as the primal problem and obtain the primal problem (5) as the dual of the dual. We do this because the dual of C(Ω, R k×d ) is known, while the dual of M(Ω, R k×d ) is difficult to characterize.
Consider the problem maximize which is equivalent to (9). Define Rewrite (10) as and consider its Fenchel-Rockafellar dual , u satisfies the zero-flux boundary condition.
We now state the rigorous form of the primal problem (5) minimize The point φ = 0 satisfies the assumption of Theorem 5. Furthermore, it is easy to verify that the optimal value of the dual problem (6) is bounded. This implies strong duality, (11) is feasible, and (11) has a solution.

Algorithmic preliminaries
Consider the Lagrangian for the vector optimal transport problems (5) and its dual (6) which is convex with respect to u and w and concave with respect to φ. Finding a saddle point of (12) is equivalent to solving (5) and (6), when the primal problem (5) has a solution, the dual problem (6) has a solution, and the optimal values of (5) and (6) [48] for further discussion on this point.
To solve the optimal transport problems, we discretize the continuous problems and apply PDHG method to solve the discretized convex-concave saddle point problem.

PDHG method
Consider the convex-concave saddle function where f , g, and h are (closed and proper) convex functions and x ∈ R n , y ∈ R m , z ∈ R l , A ∈ R l×n , and B ∈ R l×m . Note L is convex in x and y and concave in z. Assume L has a saddle point and step sizes µ, ν, τ > 0 satisfy 1 > τ µλ max (A T A) + τ νλ max (B T B).
Write · 2 for the standard Euclidean norm. Then the method converges to a saddle point. This method is called the the Primal-Dual Hybrid Gradient (PDHG) method or the (preconditioned) Chambolle-Pock method [23,10,46]. PDHG can be interpreted as a proximal point method under a certain metric [30]. The quantity is the fixed-point residual of the non-expansive mapping defined by the proximal point method. Therefore R k = 0 if and only if (x k , y k , z k ) is a saddle point of L, and R k decreases monotonically to 0, cf., review paper [49]. We can use R k as a measure of progress and as a termination criterion.

Shrink operators
As the subproblems of PDHG (13) are optimization problems themselves, PDHG is most effective when these subproblems have closed-form solutions. The problem definitions of scalar, vector, and matrix-OMT involve norms. For some, but not all, choices of norms, the "shrink" operators shrink(x 0 ; µ) = argmin for u ∈ R k×d and The shrink operators of these norms have closed-form solutions. For the matrix-OMT, we focus on norms for U ∈ H d and · 2 , · 1 , and · 1,nuc for W ∈ S , which are defined likewise. The nuclear norm · nuc is the sum of the singular values. The shrink operators of these norms have closed-form solutions. We provide further information and details on shrink operators in the appendix.

Algorithms
We now present simple and parallelizable algorithms for the OMT problems. These algorithms are, in particular, very well-suited for GPU computing. In Section 7 and 8 we deal with discretized optimization variables that approximate solutions to the continuous problems. For simplicity of notation, we use the same symbol to denote the discretizations and their continuous counterparts. Whether we are referring to the continuous variable or its discretization should be clear from context.
As mentioned in Section 6, these methods are the PDHG method applied to discretizations of the continuous problems. In the implementation, it is important to get the discretization at the boundary correct in order to respect the zero-flux boundary conditions. For interested readers, the details are provided in the appendix.
Instead of detailing the somewhat repetitive derivations of the algorithms in full, we simply show the key steps and arguments for the u update of vector-OMT. The other steps follow from similar logic.
When we discretize the primal and dual vector-OMT problems and apply PDHG to the discretized Lagrangian form of (12), we get Since the minimization splits over the i, j indices, we write At the boundary, these manipulations need special care. When we incorporate ghost cells in our discretization, these seemingly cavalier manipulations are also correct on the boundary. We further explain the ghost cells and discretization in the appendix.

Scalar-OMT algorithm
The scalar-OMT algorithm can be viewed as a special case of vector-OMT or matrix-OMT algorithms. This scalar-OMT algorithm was presented in [34], but we restate it here for completeness.

Matrix-OMT algorithm
Write shrink u and shrink w for the shrink operators with respect to · u and · w . The matrix-OMT algorithm is as follows:

Parallelization
For the vector-OMT, the computation for the u, w, and φ updates splits over the indices (i, j), i.e., the computation splits pixel-by-pixel. Furthermore, the u and w updates can be done in parallel. Parallel processors handing jobs split over (i, j) must be synchronized before and after the φ update. The same is true for the scalar-OMT and matrix-OMT. This highly parallel and regular algorithmic structure makes the proposed algorithms very well suited for CUDA GPUs. We demonstrate this through our experiments.

Examples
In this section, we provide example applications of vector and matrix-OMT with numerical tests to demonstrate the effectiveness of our algorithms. As mentioned in the introduction, potential applications of vector and matrix-OMT are broad. Here, we discuss two of the simplest applications.
We implemented the algorithm on C++ CUDA and ran it on a Nvidia GPU. For convenience, we MEXed this code, i.e., the code is made available as a Matlab function. For scientific reproducibility, we release this code.

Color images
Color images in RGB format is one of the more immediate examples of vector-valued densities. At each spatial position of a 2D color image, the color is represented as a combination of the three basic colors red (R), green (G), and blue (B). We allow any basic color to change to another basic color with cost c 1 for R to G, c 2 for R to B, and c 3 for G to B. So the graph G as described in Section 3.1 has 3 nodes and 3 edges.
Consider the two color images on the domain Ω = [0, 1] × [0, 1] with 256 × 256 discretization shown in Figures 2a and 2b. The initial and target densities λ 0 and λ 1 both have three disks at the same location, but with different colors. The optimal flux depends on the choice of norms. Figures 2c and 2d, show fluxes optimal with respect to different norms.
Whether it is optimal to spatially transport the colors or to change the colors depends on the parameters c 1 , c 2 , c 3 , α as well as the norms · u and · w . With the parameters of Figure 3a it is optimal to spatially transport the colors, while with the parameters of Figure 3b it is optimal to change the colors.
Finally, we test our algorithm on the setup of Figure 2c with grid sizes 32 × 32, 64 × 64, 128 × 128, and 256 × 256. Table 1 shows the number of iterations and runtime tested on a Nvidia Titan Xp GPU required to achieve a 10 −3 precision, measured as the ratio between the duality gap and primal objective value.

Diffusion tensor imaging
Diffusion tensor imaging (DTI) is a technique used in magnetic resonance imaging. DTI captures orientations and intensities of brain fibers at each spatial position as ellipsoids and gives us a matrix-valued density. Therefore, the metric defined by the matrix-OMT provides a natural way to compare the differences between brain diffusion images. In Figure 4 we visualize diffusion tensor images by using colors to indicate different orientations of the tensors at each voxel. In this paper, we present simple 2D examples as a proof of concept and leave actual 3D imaging for a topic of future work.   Table 2: Distances between the three images Λ 0 , Λ 1 and Λ 2 .
Consider three synthetic matrix-valued densities Λ 0 , Λ 1 , and Λ 2 in Figure 5. The densities Λ 0 and Λ 1 have mass at the same spatial location, but the ellipsoids have different shapes. The densities Λ 0 and Λ 2 have the same ellipsoids at different spatial locations.
Again, we test our algorithm on the setup of Figure 5 with α = 1 on grid sizes 32 × 32, 64 × 64, 128 × 128, and 256 × 256. Table 3 shows the number of iterations and runtime tested on a Nvidia Titan Xp GPU required to achieve a precision of 10 −3 , measured as the ratio between the duality gap and primal objective value.

Conclusions
In this paper, we studied the extensions of Wasserstein-1 optimal transport to vector and matrix-valued densities. This extension, as a tool of applied mathematics, is interesting if the mathematics is sound, if the  numerical methods are good, and if the applications are interesting. In this paper, we investigated all three concerns. From a practical viewpoint, that we can solve vector and matrix-OMT problems of realistic sizes in modest time with GPU computing is the most valuable observation. Applying our algorithms to tackle real world problems in signal/imaging processing, medical imaging, and machine learning would be interesting directions of future research.
Another interesting direction of study is quadratic regularization. In general, the solutions to vector and matrix-OMT problems are not unique. However, the regularized version of (5) is strictly convex and therefore has a unique solution. A similar regularization can be done for matrix-OMT. As discussed in [34,50], this form of regularization is particularly useful as a slight modification to the proposed method solves the regularized problem.

Acknowledgments
We would like to thank Wilfrid Gangbo for many fruitful and inspirational discussions on the related topics. The Titan Xp's used for this research were donated by the NVIDIA Corporation. This work was funded by ONR grants N000141410683, N000141210838, DOE grant de-sc00183838 and a startup funding from Iowa State University.

A Discretization
Here, we discribe the discretization for the vect-OMT problem (5) and (6). The discretization for the matrix-OMT problem is similar. We consider the 2D square domain, but (with more complicated notation) our approach immediately generalizes to more general domains. Again, we use the same symbol for the discretization and its continuous counterpart. Consider a n × n discretization of Ω with finite difference ∆x in both x and y directions. Write the x and y coordinates of the points as x 1 , . . . , x n and y 1 , . . . , y n . So we are approximating the domain Ω with {x 1 , . . . , x n } × {y 1 , . . . , y n }. Write C(x, y) for the ∆x × ∆x cube centered at (x, y), i.e., We use a finite volume approximation for λ 0 , λ 1 , w and φ. Specifically, we write λ 0 ∈ R n×n×k with λ 0 ij ≈ C(xi,yj ) λ 0 (x, y) dxdy, for i, j = 1, . . . , n. The discretizations λ 1 , φ ∈ R n×n×k and w ∈ R n×n× are defined the same way. Write u = ( u x , u y ) for both the continuous variables and their discretizations. To be clear, the subscripts of u x and u y do not denote differentiation. We use the discretization u x ∈ R (n−1)×n×k and u y ∈ R n×(n−1)×k . For i = 1, . . . , n − 1 and j = 1, . . . , n u x,ij ≈ C(xi+∆x/2,yj ) u x (x, y) dxdy, and for i = 1, . . . , n and j = 1, . . . , n − 1 u y,ij ≈ C(xi,yj +∆x/2) u y (x, y) dxdy.
In defining u x and u y , the center points are placed between the n × n grid points to make the finite difference operator symmetric.
So (∇ φ) x ∈ R (n−1)×n×k and (∇ φ) y ∈ R n×(n−1)×k , and the ∇ x is the transpose (as a matrix) of −div x . Ghost cells are convenient for both describing and implementing the method. This approach is similar to that of [34]. We redefine the variable u = ( u x , u y ) so that u x,ij = u x,ij for i < n 0 for i = n u y,ij = u y,ij for j < n 0 for j = n , for i, j = 1, . . . , n, and u x , u y ∈ R n×n×k . We also redefine ∇ φ = ((∇ φ) x , (∇ φ) y ) so that for i, j = 1, . . . , n, and (∇ φ) x , (∇ φ) y ∈ R n×n×k . With some abuse of notation, we write w ij w .

B Shrink operators with closed-form solutions
All shrink operators we consider in this paper can be built from these shrink operators. These ideas are well-known to the compressed sensing and proximal methods community [44]. However, there is a subtlety we must address when applying the shrink operators to matrix-OMT: the U update is defined as the minimization over H, not all of C, and the W update is defined as the minimization over S, not all of C. Fortunately, this is not a problem when we use unitarily invariant norms (such as the nuclear norm) thanks to the following lemma. Lemma 1. Let M ∈ H and N ∈ S. Let · be a unitarily invariant norm, i.e., U AV = A for any A ∈ C, and U, V ∈ C unitary. Then is Skew-Hermitian.
Proof. The minimum over H ⊂ C is the same as the minimum over all of C, if the minimum over all of C is in H. Likewise, the minimum over S ⊂ C is the same as the minimum over all of C, if the minimum over all of C is in S.
Write M = U ΣU T for M 's eigenvalue decomposition, which is also its singular value decomposition. Define σ + = argmin s {f (s) + (1/2) s − σ(M ) 2 2 }, which implies 0 ∈ ∂f (σ + ) + σ + − σ(M ). With this, we can verify that 0 ∈ U diag(∂f (σ + ))U T + U diag(σ + )U T − U diag(σ(M ))U T , i.e., U diag(σ + )U T satisfies the optimality conditions of the optimization problem that defines M + . So M + = U diag(σ + )U T , which is Hermitian. Likewise, write N = U ΛU T for N 's eigenvalue decomposition. N has orthonormal eigenvectors and its eigenvalues are purely imaginary. We separate out the magnitude and phase of Λ and write N = U diag(σ(N ))P U T . More precisely, σ(N ) = | diag(Λ)| and P is a diagonal matrix with diagonal components ±i. Then the SVD of N is N = U diag(σ(N ))V T where P U T = V T . With the same argument as before, we conclude that N + = U diag(σ + )V T for some σ + ∈ R n . So N + = U Λ + U T where Λ + = diag(σ + )P is diagonal and purely imaginary, and we conclude N + is skew-Hermitian.
The norm · 1 as described in Section 6.2, is not unitarily invariant. However, its shrink operator acts element-by-element, so it is easy to arrive at a conclusion similar to that of Lemma 1.


PanJen: An R package for Ranking Transformations in a Linear Regression

PanJen is an R-package for ranking transformations in linear regressions. It provides users with the ability to explore the relationship between a dependent variable and its independent variables. The package offers an easy and data-driven way to choose a functional form in multiple linear regression models by comparing a range of parametric transformations. The parametric functional forms are benchmarked against each other and a non-parametric transformation. The package allows users to generate plots that show the relation between a covariate and the dependent variable. Furthermore, PanJen will enable users to specify speciﬁc functional transformations, driven by a priori and theory-based hypotheses. The package supplies both model ﬁts and plots that allow users to make informed choices on the functional forms in their regression. We show that the ranking in PanJen outperforms the Box-Tidwell transformation, especially in the presence of inefﬁciency, heteroscedasticity or endogeneity.

Introduction
A model that fits data well but is unrelated to theory can only describe correlations. In contrast, a model with both a good fit and a theoretically sound foundation can give insights into hypotheses on causality. The functional form in a regression model describes the relationship between a dependent variable and its covariates. There are numerous examples of researchers who have neglected functional form relationships and applied the default linear relationship between variables in their regression models (Box, 1976;Breiman and others, 2001;Berk, 2004;Angrist and Pischke, 2010). From a superficial point of view, these models may provide efficient parameter estimates with narrow standard errors, high t-values, and significance. A strong non-linear relationship between a dependent variable and a covariate will often offer reasonable test statistics with a default linear functional form specification. However, positive test statistics are not the same as proof of a linear relationship. Even more important is the fact that a misspecified functional form can lead to an incorrect interpretation and prediction of the relationship between the dependent variable and a given covariate. The specification should be driven by theory with an a priori hypothesis of the relationship between the dependent and covariates.
PanJen was developed over several years of applied research on property value models. Here, the sales price is estimated as a function of its characteristics, such as the size of the living space, the number of rooms, and access to shopping. However, the package is applicable to most cases in which the relationship between a continuous dependent variable and its covariates is explored. In the applied econometrics literature, this task has commonly been solved using power transformations in initial analyses (Palmquist, 2006). Two examples of power-transformations are Box-Cox and Box-Tidwell (Box and Tidwell, 1962;Box, 1976). The power-transformations are easy to use; the ability of these two power transformations to detect functional forms has been studied extensively in the academic literature, e.g., Kowalski and Colwell (1986); Brennan et al. (1984); Clark (1984), and they are still used in applied studies (Cohen et al., 2013;Farooq et al., 2010;Link, 2014;Joshi et al., 2017;Benson et al., 1998;Troy and Grove, 2008).
The popularity of these power-transformations is surprising given that their shortcomings are well described in the literature (Levin et al., 1993;Wooldridge, 1992a). While power transformations perform well in many circumstances, they do not perform well in the presence of omitted variables, inefficiency, heteroscedasticity and endogeneity. Furthermore, a transformation can be challenging to interpret, does not necessarily relate back to a theory-driven hypothesis and does not detect whether the relationship changes across the distribution of the dependent variable.
Another approach to the functional form issue is to abandon the parametric model and approach the challenge from a non-or semi-parametric angle. In the academic literature, a number of alternatives have been proposed and used, such as non-parametric or semi-parametric methods (Anglin and Gençay, 1996;Gençay, 1996;Clapp and Giaccotto, 2002;Bin, 2004;Geniaux and Napoléone, 2008). Non-or semi-parametric models provide data-driven approaches to establish the relationship between a dependent variable and covariates. A non-parametric model can be attractive, because the functional form is revealed by the data instead of being be predefined by the researcher. However, the gained flexibility of non-parametric analysis comes at the cost of more difficult interpretation of the estimates, which is perhaps why parametric models are often used in applied work.
In non-parametric models, the relationship is fitted to the sample to the extent that the estimated relationship is at risk of being over-fitted. In other words, the estimated effect captures random error or noise in combination with the underlying relationship in the population (Wood, 2006). Additionally, a central critique concerning this approach is that the results from a non-parametric model are difficult to generalize or extend outside of the sample (McMillen and Redfearn, 2010). Even so, a non-parametric model holds great potential in exploratory analysis.
We are certainly not the first researchers to consider the possibility of utilizing the apparent advantages of the non-parametric modelling to explore functional form relationships in parametric modelling. The literature on using a non-parametric model to test different parametric specifications is large (González-Manteiga and Crujeiras, 2013). The primary approach in the existing literature has been to test a parametric version of a model against a non-parametric version (Wooldridge, 1992b;Horowitz and Härdle, 1994;Zheng, 1996;Li et al., 2016). However, to the best of our knowledge, none of these tests have been widely adopted in the empirical literature. PanJenis developed to allow applied researchers to utilize non-parametric estimation to identify a better parametric functional form. The package offers a test based on well-established measures and is provided on a well-established software platform. We believe that very few empirical researchers know about the existing tests, and the few that do perceive them as too complicated due to their non-parametric basis. PanJen offers the user a transformation-ranking based on the parametric transformation that captures most of the variance of the dependent variable. The ranking includes a non-parametric specification that can detect if the relationship between the dependent variable and covariate is non-stationary. In contrast to existing tests, semi-parametric transformations are only included as a benchmark rather than as an incremental part of the test. The engine in the PanJen ranking is the well-known Akaike Information Criterion (AIC ) and Bayesian Information Criterion (BIC ) measures. These model-fit measures are already widely applied in the empirical literature and should not be a hindrance for the empirical researcher. With the PanJen ranking, we hope to introduce an approach that will make applied researchers explicitly consider functional form in their parametric models.
In the next section, we briefly describe the main idea behind the ranking in the package (PanJen ranking. In the following section, we describe the workhorse behind it, the Generalized Additive Model (GAM ). In section 4, we explain how the PanJen ranking works. In section 5, we illustrate how to use the package by using a real example from our research. Section 6 offers a comparison between the PanJen ranking and the Box-Tidwell transformation. We simulated 10,000 datasets and recovered the functional form of one variable in a model with different impediments to show the merits of PanJen relative to a conventional approach. In section 6, we conclude the paper with a short discussion of when and how the researcher should use the package with an emphasis on the risk of pre-test bias.

The main idea of the PanJen ranking
PanJen is built on the idea that the choice of a functional form can be guided by model fit. In the PanJen ranking, a given number of models that vary only in the transformation of one covariate are estimated. One of these transformations is a so-called function that for now we will simply note makes this one model semi-parametric. All the estimated models are then ranked according to their BIC). The BIC provides a relative goodness-of-fit measure that accounts for the complexity of the model. More formally, the PanJen ranking estimates a model Y = β 0 + Xβ k + g(x)β l + ε where Y is the dependent variable, ε is an i.i.d. error-term, X is a vector of k of covariates, β k is the corresponding vector of parameter estimates, g(x) represents a set of functional form transformations among the set: where β l the corresponding l parameter estimates for the parametric transformations. In the last two transformations, there is no parameter estimate for g(x), because f (x) is the non-parametric smoothing and 0 leaves out the explanatory variable.
The ranked BIC-values show how each transformation performs relative to the others. The semiparametric transformation allows the user to assess how well parametric transformations perform relative to a flexible semi-parametric function. If the data generation process does not resemble any of the parametric transformations, the smoothing function will still capture the relationship. The BIC scores are supplemented by the closely related AIC. In practice, both the AIC and the BIC penalize model complexity, although the penalty term in BIC is larger than in AIC (Burnham and Anderson, 2004). The smoothing function is highly flexible, but the flexibility comes at a cost. Therefore, it is not necessarily ranked the highest since both AIC and BIC penalize the model complexity. There is no objective and transparent way to choose between the measures. The right measure depends on the user's a priori theory of the data generation process. If the users assume that one of their models perfectly fits the underlying data generation process, then the BIC is the right measure. If they instead assume that the underlying data generation process is extremely complex and none of the possible models will be able to perfectly capture it, then AIC is the right measure (Aho et al., 2014).
The PanJen ranking is supported by a plot function that graphically outlines the relationship between the dependent variable and covariate. The plot is created by predicting the dependent variable using the median for all independent variables other than the one in question. The variable in question varies across a scale from the 5th quantile to the 95th quantile of the actual distribution in the dataset. The plot shows the user how each transformation captures the relationship across the distribution of the dependent variable. If the smoothing far outperforms all parametric transformations, the reason may be that the relationship changes across the distribution and the proposed simple parametric transformation does not capture the relationship between the dependent and independent variable. The plot will reveal this.

A semi-parametric model for benchmark
We estimate the parametric transformations using the Generalized Linear Model (GLM ) and the semiparametric using GAM. GAM is a special case of the Generalized Linear Model (GLM) in which it is possible to include one or more so-called smoothing functions. A smoothing function is a nonparametric way to include a continuous covariate in a parametric model and make it semi-parametric.
The GAM can be written as follows: Y i is the dependent variable of observation i. It is distributed as an exponential family distribution, e.g.the normal, the gamma or chi-square distribution. Xi is a matrix of covariates that are parametrically related to the dependent variable. β is the corresponding vector of the parameter estimate, and f i is a smoothing function of covariate x 1i .
The GAM provides a flexible specification of a covariate by only specifying it as a smoothed function. By entering a variable with a smoothing function, the researcher does not specify a functional form, but instead lets the data speak. The smoothing function comprises the sum of k thin plate regression spline bases b h (•) multiplied by their coefficients. It is estimated as follows: . The non-parametric component of the model is fitted with a penalty on wiggliness (how flexible the smoothing is). The penalty, θ, is determined from the data using generalized cross-validation or related techniques. The penalty directly enters the objective function through an additional term capturing wiggliness in the smoothing function, i.e., Here, Y is the fitted dependent variable, and the second derivatives of the smoothing function describe its wiggliness. We estimate the GAM using the mgcv R-package mgcv (Wood, 2017). For a thorough introduction to GAM, please see (Wood, 2006).

Using the package
We illustrate the use of PanJen using a hedonic house pricing model. The central idea is that the sale price of a home is a function of its characteristics, understood as both the characteristics of the home itself and its surroundings. The latter poses a problem in the empirical application of the hedonic method because observations can be correlated through space. A very flexible solution to this problem is to use the GAM framework to smooth over the x-y coordinates, thus allowing one to non-parametrically control for spatial correlations. von Graevenitz and Panduro (2015) illustrated the relationship between smoothing over space and classic spatial econometrics with weight-matrices and fixed spatial effects. They also showed that smoothing is a better alternative when the researcher does not know the underlying spatial data generation process. For recent applications, please see Rajapaksa et al. (2017) or Schäfer et al. (2017).
In our example here, we solely focus on the structural characteristics. These characteristics are measured by a range of variables. The researcher does not a priori know how the characteristics of the house are related to its price. For example, we expect the price to increase with the size of the home, but we do not know if that relationship is linear. It could be that going from 2 to 3 bedrooms is different than going from 7 to 8 bedrooms, i.e., we want to know if we should take account of marginally increasing or decreasing price-relationships. PanJen was developed to answer this type of question by finding the functional form relationship between the home price and different home characteristics.

An example: the implicit price for living area
Names Description lprice log transformed price in 1000 EUR area living area in square meters age build year bathrooms number of bathrooms lake_SLD distance to nearest lake in meters highways distance to nearest highway in meters big_roads distance to nearest large road in meters railways distance to nearest railway in meters nature_SLD distance to nearest nature area in meters  (Lundhede et al., 2013). We have 9 continuous and 7 dummy variables for quality at our disposal. In addition, the dataset includes 3 year dummies to control for price trends. The variables are listed in Tables 1 and 2  > library(PanJen) > data("hvidovre") Then, we set up a formula-object. We log-transform the prices because this introduces flexibility and is the convention within the hedonic literature (Diewert, 2003). It is possible to test different transformations by simply transforming the variable or test different link-functions by leaving variable empty in fform(). Ten of the variables are dummy variables where transformations are irrelevant. We include only these in the first regression: > formBase<-formula(lprice~brick+roof_tile+roof_cemen + + rebuild70+rebuild80+rebuild90+rebuild00+y7+y8+y9) > summary(gam(formBase, method="GCV.Cp",data=hvidovre)) Family: gaussian Link function: identity Formula: lprice~brick + roof_tile + roof_cemen + rebuild70 + rebuild80 + rebuild90 + rebuild00 + y7 + y8 + y9 Parametric coefficients:  This initial model explains nearly 15% of the variation in price. The next characteristic we want to control for is the size of the home. In the dataset, the living area in square metres is stored under 'area'.
We start out by using the default transformations supplied by the PanJen function fform(). This function ranks the fit of nine predefined transformations and a smoothing. The mandatory inputs are the name of the dataset, the model formula and the new variable we wish to test using the PanJen ranking: > PanJenArea<-fform(hvidovre,"area",formBase,distribution=Gamma (link=log) Wood (2006) for an elaboration" The results are ranked according to their BIC. Strictly according to this ranking, we should logtransform the area. This implies approximately that a % change in living area results in a % change in price. Given the respondent variable has been log-transformed previous to the model-fitting, any interpretation at the the original scale should be done with care , see e.g. (Barrera-Gómez et al., 2015).
The differences in the score for the four lowest BIC are small, and it might be a matter of differences in the tails of the distribution. This can be checked by plotting the predicted price against the area. The function plotff() generates a plot with the predicted price against the area from the 5th to the 95th percentiles with all other covariates variables at their median value:

> plotff(PanJenArea)
The black line is the smoothing function. The log-squared and the linear specification closely follow this line. In conclusion, the implicit price for the living area is positive and slightly marginally declining. You can specify your own transformations using choose.fform(). In the following, we test three transformations: 'area', 'log(area)' and 'area 2 '. We start out by defining a list of transformations:

A changing relationship
The first transformation was rather straightforward because the relation between area and price is stable across the price-distribution. The age of the home is a characteristic with which this is not the case. The reason is that the building changes over the years. As the home gets older, there is not a direct link between the age and state of the home. Where a newly built home is likely to be built with modern standards and tastes in mind, an old house can instead be charming and authentic. At the same time, in general, houses built during periods of building booms, which in Denmark were in the 1960s, are of lesser quality than those built in the 1950s or 1970s. We can show this by running fform() and plotting the results.
The plot shows the relationship between age and price from the 5th to 95th percentile of the age Figure 3: The relation between price and age distribution. Given the plot, it is difficult to think of a parametric relationship that will capture this relationship. If the age of the home is somehow related to the research question, the best solution might be to use the smoothing function. If age is nothing more than a control variable, one could perhaps resort to interval dummies similar to the year dummies in the model. As a part of the final test of the model, it would be worthwhile to test to what degree the variable of interest is robust to the way age enters the pricing function. In our setting, what should be noted is that the complexity of the relation between age and price would have gone unnoticed if we had compared only the parametric transformations.

Interacting with other packages:
When you run choose.fform() or fform(), all generated models and datasets are stored in a new list of list -object. Within the list 'models', all estimated models are stored as 'gam ', 'glm' and 'lm'-objects. This means that all objects used to create the fform output are easily available. The plotting function in PanJen is simple, but perhaps not enough when the researcher needs to produce plots for a third party. Here, we show how to use this to make a plot using base R, but we could just as well have used the tlm if we wanted a more detailed plot (Barrera-Gómez and Basagaña, 2017).

Monte Carlo simulations
We tested the performance of the PanJen ranking against the well-known Box-Tidwell transformation (Box and Tidwell, 1962). We choose this as a benchmark since it is the most structured choice applied in existing empirical work. Power transformations such as Box-Cox and Box-Tidwell were suggested in the 1960s by Box and Cox (1964) and Box and Tidwell (1962). The Box-Tidwell transformation identifies the transformation that minimizes non-normality in the error term and linearizes the relationship between the dependent variable and the covariate using a maximum likelihood function. Thus, the researcher can use the test to find the power-transformation with the highest likelihood. This section presents the results from nine Monte Carlo simulations in which the performance of the Box-Tidwell and PanJen is tested.
The simulations are centred on a base model: where x 3 is the variable of interest and x 1 and x 2 are two other covariates. The functional relationship between Y and f (x 3 ) was then tested using PanJen ranking and Box-Tidwell transformations. Table 3 summarizes the results. The fourth and fifth columns show the share of times each method reported the true functional form. In the Box-Tidwell case, the transformation parameter was allowed to vary by up to 0.2 from the correct specification.

Simulation Simulation description
PanJen BoxTidwell  Each of the nine simulations tested the robustness of the methods in relation to different wellknown econometric methods. Overall, the PanJen ranking performed acceptably. The method pointed to the correct functional form in 97 to 100% of the cases. The Box-Tidwell transformation performed just as well when the dataset was "well-behaved." It was already well-established in the literature that the method is sensitive to endogeneity, inefficient model estimates, heteroscedasticity and endogeneity, and this is also what we find in our study. In conclusion, PanJen Ranking performs better or just as well as Box-Tidwell.

Conclusion
In this paper, we present the PanJen package. We provide a simple and intuitive description of the PanJen ranking. Based on a house price dataset, we show how the functions in the package can be applied to determine the relationship between a dependent variable and its covariates. Furthermore, we compare the PanJen ranking method to the Box-Tidwell transformation and show that the PanJen ranking performs just as well as or better than the Box-Tidwell transformations. The Panjen ranking outperforms Box-Tidwell in situations where the model suffers from inefficiency, heteroscedasticity or endogeneity. In some circumstances, the theory provides little or no guidance on the functional relationship between the dependent and covariates in multiple regression models. In such circumstances, PanJen can support users in their decision on the functional form of the covariates. If the functional form relationship is more complex than a simple parametric transformation, we suggest considering a semi-or non-parametric model. The package has deliberately been restricted to test one covariate at a time without a silent output option. We want to deter the user from looping over every explanatory variable in search of a fit using the PanJen ranking, because this increases the pre-test bias. However, we also recognize that exploratory analysis is part of any empirical application of statistical modelling. Learning is a sequential process, and in many circumstances, we have not properly thought out an a priori hypothesis on which to base our models (Wallace, 1977). People perform exploratory model estimations and in many cases under-report their approach. Even so, pre-test bias is not a problem caused or exacerbated by PanJen Ranking. Regardless of how a researcher performs multiple model estimations, the risk of pre-test bias can be reduced by adopting a sampling approach. The sampling approach can be implemented by dividing data into training and test datasets, where the explorative analysis is conducted on the first. It is our hope that people will use PanJen to improve their models by specifying relationships that more accurately fit their data. In doing so, users should consider PanJen ranking as a guide and not as a substitute for a priori hypothesis.