Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
96,789 characters · 33 sections · 47 citation commands
Machine Learning based Framework for Robust Price-Sensitivity Estimation with Application to Airline Pricing
\abstract{Problem definition: We consider the problem of dynamic pricing of a product in the presence of feature-dependent price sensitivity. Developing practical algorithms that can estimate price elasticities robustly, especially when information about no purchases (losses) is not available, to drive such automated pricing systems is a challenge faced by many industries. Methodology: Based on the Poisson semi-parametric approach, we construct a flexible yet interpretable demand model where the price related part is parametric while the remaining (nuisance) part of the model is non-parametric and can be modeled via sophisticated machine learning (ML) techniques. The estimation of price-sensitivity parameters of this model via direct one-stage regression techniques may lead to biased estimates due to regularization. To address this concern, we propose a two-stage estimation methodology which makes the estimation of the price-sensitivity parameters robust to biases in the estimators of the nuisance parameters of the model. In the first-stage we construct estimators of observed purchases and prices given the feature vector using sophisticated ML estimators such as deep neural networks. Utilizing the estimators from the first-stage, in the second-stage we leverage a Bayesian dynamic generalized linear model to estimate the price-sensitivity parameters. Results: We test the performance of the proposed estimation schemes on simulated and real sales transaction data from the Airline industry. Our numerical studies demonstrate that our proposed two-stage approach reduces the estimation error in price-sensitivity parameters from 25% to 4% in realistic simulation settings. \textbf{\textit{Managerial implications}}: The two-stage estimation techniques proposed in this work allows practitioners to leverage modern ML techniques to robustly estimate price-sensitivities while still maintaining interpretability and allowing ease of validation of its various constituent parts.}
Keywords: Dynamic Pricing, Price Elasticity, Causal Inference, Machine Learning
Many sellers are interested in improving their pricing decisions by dynamically adjusting prices of their products based on unique product features and other relevant information available at the time of request e.g., competing product prices, market indices etc. Automated dynamic pricing systems that can enable such functionality typically require estimates of customers' price-sensitivity at this granular level using historical sales transaction data. Although modern IT systems have made it possible to collect and use rich request-specific information from a shopping session, the estimation of price-sensitivity from such feature-rich historical observational data still remains a challenging task in many industries hua2021kdd, ye2018kdd. This is especially true when the firm is only able to observe the total quantity sold for a product at a certain price in given time period but is unable to observe customer requests that did not result in a sale. Data about unsuccessful sale events at an individual request level is commonly referred to as loss or no-purchase information and if available can be utilized to estimate price-sensitivity by directly modeling purchase probabilities as a function of price and other relevant features via methods like parameterized logistic regression. Situations where loss information is not available are common in industries like Airlines, Car Rentals, and Hotels where a large fraction of sales occur via third-party websites or Online Travel Agents (OTA) such as Expedia, KAYAK, Google, Orbitz etc. Typically, a firm selling its products via these OTAs does not have direct visibility into loss information corresponding to a large fraction of individual requests on these platforms. In the absence of such loss information, one may estimate price-sensitivity by learning price-demand relationship (demand response modeling) using total sales for a product in given time period e.g., daily sales.
In this work we focus on learning such price-demand relationship given historical observations. This is akin to an observational study setup where the purchases are observed under the historical pricing policy deployed by the firm. A na\"ive approach can be to fit a supervised model with desired features (including price-related) and demand as the response variable. But this approach has several shortcomings. One of the shortcomings is related to the choice of the model. On one hand, parametric models that are popular in the literature provide an explainable and easy to use model with interpretability of their parameters, but are limited in terms of their flexibility and have strong assumptions about the generating process. On the other hand, modern machine learning (ML) approaches such as neural networks and gradient boosting, despite their high predictive power, do not easily lend themselves to constructing an interpretable framework for the inference task related to the price sensitivity estimation.
Another issue of directly fitting a supervised model is the presence of confounders, which are features that affect both demand and price. It is well known that omission of such features from the regression model may lead to biased estimates of price sensitivity parameters when estimated from observational data. In addition, although our main goal in this paper is not demand prediction, it is worth mentioning that a complex predictive model directly fitted on demand using price and additional features cannot generally make reasonable demand predictions for prices outside the range of observed values or prices different from the pricing policy underlying the observed data, as those regions are out of the training data distribution zou2020counterfactual,shalit2017estimating. Explainable and simpler parametric models based on micro-economic theory and domain specific knowledge may be less susceptible to this challenge compared with more complex approaches such as deep neural networks.
Our approach for demand response modeling is thus based on a hybrid framework where we utilize modern ML-based non-parametric techniques such as deep neural networks to make an observational baseline prediction of both price and demand, and then combine it with a parametric approach to estimate the causal relationship between deviations of price and demand from this baseline. To this end, we build upon recent advances in the area of causal inference and semi-parametric regression in econometrics Athey2019 and employ related ideas to construct robust price-sensitivity estimators using Poisson semi-parametric modeling. Moreover, by utilizing variational Bayesian approach for updating the model parameters in an online fashion, this approach is compute and storage efficient from a practical implementation perspective.
In the airline dynamic pricing application motivating this work, a product being priced is typically an itinerary defined by specific route (Origin, Destination and Path), trip type (round trip, one way, multi-city etc), date and time of departure. Moreover, airlines typically differentiate products by compartment (economy, economy-premier, business, first etc.) and associated privileges and restriction e.g., change-related fees, upgrade eligibility, baggage allowance etc. Additional trip-related features are also used to differentiate customer price-sensitivity, these may include duration of stay at the destination, advance purchase restrictions, departure day of week, and number of layovers.
Airline demand is significantly impacted by yearly seasonal patterns; for example often there is more demand in summer months due to school holidays. Moreover, special events (HSE) such as conferences in the origin or destination city also have a major impact on demand. Demand patterns also vary over the booking horizon of an itinerary. Most airlines allow the booking process to start around one year before the departure date. The pattern of demand over this booking horizon can show significant variation due to the fact that leisure customers, who represent a large fraction of demand, are highly price-sensitive, and book their trips early while business travelers with higher willingness-to-pay (WTP) may book their trip much later. Therefore, causal inference of price-sensitivity parameters via price-demand relationships requires inclusion of features to control for yearly seasonality, booking days prior and HSEs in the model.
The competitive landscape that the airline is operating in also plays an important role in determining its demand. On most routes there are multiple airlines offering services. The set of competing products that may impact an itinerary's demand are not only the similar (and relevant) itineraries offered by a competing airline but also other products for the same origin-destination pair offered by airline itself. Moreover, quality indicators such as on-time performance or rating of in-flight services could also be correlated to price-sensitivity.
Unlike many scenarios related to e-commerce setting where the inventory costs for the product being sold are minimal due to supply being relatively elastic, airlines have to deal with fixed supply of seat inventory on flights in their network. Any seat not sold by the departure date brings no revenue and gets wasted i.e., the inventory is perishable. To deal with the fixed and perishable inventory, airlines employ sophisticated Revenue Management (RM) systems. Prices generated by an airline for an itinerary are the product of complex interactions between the RM system output and further refinements done via pricing rules, also known as fare rules in the airline industry. Airline RM systems work at a higher (aggregated) level than the level at which an itinerary is priced. For example, regardless of whether an RM system performs optimization at the leg or the origin-destination (OD) level, round-trip related dimensions e.g., duration of stay, are not explicitly considered in the RM system optimization process. RM systems determine the expected marginal opportunity costs involved in fulfilling the itinerary request by considering product prices at an aggregated level (fare classes), forecasts of future demand, remaining seat inventories and time until departure. These marginal opportunity costs, also known as bid prices in the airline industry, are used to determine availability of fare classes at the OD level. Further refinement to get to the level of itinerary price is handled via a complex set of fare rules tied to these fare classes Vinod2021. Often these fare rules are not fully data-driven based on an empirical understanding of customers' price-sensitivities, but rely to a large extent on experience of analysts and managers in the airline's pricing departments.
In this work, we focus on developing more data-driven approaches to refine prices given the aggregate level marginal opportunity costs (bid prices) computed by the RM system. In particular, we assume that the firm has an RM system to perform finite-horizon network level optimization at an aggregate level and therefore it can compute the marginal opportunity costs, $c_{it}$, associated with fulfilling an itinerary request $i$ that it receives at time $t$. Given this cost and an $m$-dimensional feature-vector, $\mathbf{X}_{it}$, associated with this itinerary request, the firm is interested in further refining the price to charge for this request. Features associated with the itinerary can include information about origin, destination, number of stops, point of sale, departure time, day of week, holiday indicators etc.
To enable such dynamic pricing policies, we develop robust methods for estimating the price-sensitivity of customer population for such featurized itineraries via demand response modeling. Note that this system for further price refinement is separate from the RM system and only takes as input marginal opportunity costs generated by the RM system i.e., the problem of interest in this work is to find an optimal price for the itinerary given that the marginal opportunity cost has already been computed by the RM system. The RM and pricing systems are implicitly connected as the bookings generated and prices associated with sales transactions are fed back into the RM system forecaster at an aggregate level to forecast class level demands and fares used for network optimization. We assume that the firm takes into consideration some consistency requirements between the two systems while generating inputs for the RM system e.g., by appropriately mapping booking and price transactions to aggregated level products (classes) used in the RM system, pro-rating fares from round-trip to half-return level etc. This interaction between the RM and Pricing systems although important are beyond the scope of this work and here we solely focus on the price sensitivity estimation problem.
Given the complex RM system and pricing methods used by the airline, the price offered for an itinerary is correlated with many of the features we have discussed earlier in this section. What makes estimating price-sensitivity parameters via demand response modeling challenging is that the demand for an itinerary is also correlated with many of the same features that impact its price. For example, higher demand in summer months leads to higher prices due to the opportunity costs (or bid prices) for that time period being higher. Such features that affect both the price (treatment variable) and sales (outcome variable) are known in econometrics literature as confounders. If the confounders are not identified and included in the model i.e., they are not controlled for, a na\"ive estimation of price-sensitivity coefficients will be biased due to the so-called “omitted variable bias” Angrist2009. This bias problem cannot be solved just by using more sophisticated ML based estimators. Instead, one would need to rely on Instrumental Variable (IV) techniques or conduct randomized pricing experiments by assigning randomly chosen prices to itinerary requests Angrist2009.
Randomized pricing-based studies to estimate customers' price-sensitivity in the airline industry may also not be possible. The main reason is that such experiments could be perceived to be quite costly in the short-term because of using “non-optimal” pricing policy and fear of losing revenue. Besides, if two booking requests with the same features are offered randomly different prices, the passengers may try to get around the randomization by trying several times to get a lower price. Even worse, the randomization approach can be perceived as unfair and dishonest.
Given these challenges, industries like airlines have found it difficult to implement a robust framework for automated dynamic pricing. The current state of affairs has necessitated the use of complex rules-based pricing systems as the de facto solution to dynamic pricing in this industry. In our view, sophisticated ML techniques by themselves cannot solve this problem. However, combining these techniques with the recent advances in the area of causal inference allows for a more data-driven solution to this problem.
This work is related to the literature on dynamic pricing in the presence of contextual features, and the literature on causal inference in the presence of treatment heterogeneity with applications to price elasticity estimation.
Constructing price-demand relationships for goods has a long history in economics, specifically in macro-economic theories related to the study of supply and demand. The use of demand response modeling in the more applied area of dynamic pricing via maximizing profit or revenue functions is relatively more recent but has already lead to a substantial amount of literature. The survey paper by denBoer2015 provides an extensive overview of the dynamic pricing and learning literature. The majority of the work in this area focuses on parametric models, where the functional form of the demand function is assumed to be known. This includes more recent works on dynamic pricing in the presence of contextual features Cohen2016, Qiang2016, Javanmard2019a, ban_personalized_2021, Elmachtoub2021. There is also some work in the area of non-parametric models which do not assume a known functional form for the demand model, see for example Besbes2009, Besbes2015 and Perakis2019. However, most of these works have not focused on robust estimators for price-sensitivity from the causal inference perspective in the presence of many confounding variables, a situation common in industries like Airlines. There has also been work on semi-parametric models especially in the econometrics literature Robinson1988, Newey2004. These models assume a parametric (often linear) structure for treatment effect and a non-parametric structure for the remaining part of the model. However, most of this line of work has focused on additive noise structure for the model of interest and results for demand models using discrete distributions like Poisson or Negative Binomial has been largely lacking. Moreover, there are new challenges related to regularization induced bias in the estimation of parameters of interest if the non-parametric part is modeled using modern ML techniques chernozhukov2018double.
For estimating treatment effects e.g., price sensitivities, Instrumental Variable (IV) based two-stage estimation schemes are often used in the causal inference literature to reduce estimation bias pearl2009causality,newey2003instrumental,hartford2017deep,kmenta2010mostly. An IV is a variable that directly affects the treatment (e.g., price) but not the outcome (e.g., bookings), i.e., its effect on the outcome is only through the treatment. IVs help in estimating the causal relationship from regression models in the presence of endogenous regressors, for example due to omitted confounders. Two-stage least squares (2SLS) estimation angrist1996identification is one of the most widely known IV techniques which has a linear and homogeneous treatment effect assumption. In the first stage of 2SLS, the endogenous variable is regressed on the IVs to create a new latent variable, and in the second stage the endogenous variable is replaced by the new variable (prediction from the first stage model) in the main linear regression model of the outcome. Non-parametric IV frameworks have been developed for more general settings than the linear model newey2003instrumental,darolles2011nonparametric,chen2012estimation,hartford2017deep. In hartford2017deep, the authors train deep neural networks to model the distribution of treatment given the IVs and other features in the first stage and the outcome prediction model in the second stage. Although IV based frameworks are promising, identifying IVs itself is a challenging problem which limits their application.
With firms having access to richer data, it is possible to include most of the important confounders into the regression model. However, this results in a high-dimensional setting and may also require a more sophisticated regression model. For such settings, Double or Orthogonal machine learning (DML) is another line of causal inference methods that use a two-stage estimation process chernozhukov2018double. This approach assumes a partially linear model where the treatment effect is linear but other features are either high-dimensional (with respect to number of samples) and/or their effect on treatment and outcome is non-parametric. It combines sophisticated machine learning techniques for constructing estimators for the non-parametric part in the first-stage with a second stage estimation of the treatment effect based on Neyman-orthogonal moment equations for robust estimation of treatment effects. DML procedure can achieve a $\sqrt{n}$-consistent estimator of the treatment effect under very weak conditions. semenova2017estimation extend DML to a panel data setting with a high-dimensional treatment effect and unobserved unit heterogeneity, but limit the relationship between features and treatment/outcome to be linear. Other extensions and variations of the DML framework are also available in the literature mackey2018orthogonal,oprescu2019orthogonal. Nevertheless, DML methods and its variants have practically been mostly applicable to partially linear regression models or setups with binary treatment. nekipelov2022regularised is a recent work that proposes an Orthogonal ML approach for nonlinear semi-parametric models and is applicable to our problem setup.
Machine learning methods such as random forests and neural networks have recently been used to build methods for estimating treatment effects and counterfactual prediction. Some recent works have proposed methods for estimating non-parametric heterogeneous treatment effects based on random forests athey2019generalized,wager2018estimation. Other works have employed neural networks for learning treatment invariant representation of confounders to adjust the observational distribution shalit2017estimating, and variational autoencoders (VAE) kingma2013auto,rezende2014stochastic for inferring unobserved confounders with proxies louizos2017causal or inferring latent factors of treatment and decorrelating them with confounders zou2020counterfactual. However, all of these works focus mainly on binary (or vector of binary) treatments. In contrast, our focus in this work is on a continuous treatment variable such as price with the response being modeled using discrete distributions such as Poisson with a non-linear demand-response function to model the price-demand relationship.
In this article, we consider a setting in which a firm can measure relevant confounders for estimating price-sensitivity. However, because the firm is dealing with a large number of features with complex interactions, fully parametric models for constructing demand response functions may not suffice. To address this situation, we formulate a Poisson semi-parametric demand response model with feature-dependent price-sensitivity. This model allows us to leverage the flexibility provided by modern ML techniques, while keeping the specification of the effect of price on demand interpretable. Pricing methodologies can have a large impact on the firm's revenue and for such critical tasks, interpretability vis-a-vis price sensitivity parameters is important as it allows a firm to qualitatively validate the inferred parameters easily by considering the individual effect of each feature. For example, parameters associated with day of week can be used to analyze relative impact of each day on price-sensitivity to see if weekdays have higher willingness-to-pay as compared to weekends on markets with predominant business traffic. Similarly, parameters associated with time of day can be checked to see if morning and evening flights have relatively higher willingness-to-pay as compared to afternoon flights on weekdays etc. We propose two approaches to construct estimators of price sensitivity parameters associated with this demand model.
The modeling and estimation aspects related to the Poisson semi-parametric model are not straightforward. To address that, we first develop a direct approach to implement this model using deep neural networks. The estimates for price-sensitivities can be directly obtained from the weights of one of the layers of the trained neural network. While this methodology is relatively easy to implement, it may be challenging to validate, thereby making it less feasible from a practical application perspective. Moreover, the estimates of price-sensitivities may be susceptible to bias since we simultaneously estimate the price-sensitivity parameters with other unknown parameters of the model and use regularization techniques to prevent overfitting.
To address some of the issues of the direct estimation approach described above, our main contribution involves development of a robust two-stage approach for estimating price-sensitivities by extending the ideas related to Double or Orthogonal machine learning chernozhukov2018double to a Poisson semi-parametric setting. In the first stage of the method, we leverage the predictive power of modern machine learning methodologies to build predictive models for price and demand conditional on other features. In the second stage, we utilize interpretable parametric models and employ dynamic Bayesian models with an efficient sequential learning scheme based on variational Bayes ideas and Laplace approximation. The Bayesian model is naturally equipped with uncertainty estimation, allows for encoding prior information, and exhibits robustness in our experiments. In addition, it is suitable for streaming data and works well in the presence of time-variation of parameters that may exist in practice.
We also discuss the related price optimization problem and how the uncertainty estimates related to price-sensitivities can be used for evaluating policies for price experimentation based on Reinforcement Learning (RL) techniques such as Upper Confidence Bound (UCB) etc. Sutton2018.
To compare the performance of the proposed methodologies, we conduct a simulation study and show that the two-stage construction leads to robust and more accurate estimators of price-sensitivity parameters when compared to the direct one-stage approach. We also conduct a numerical study using real data to qualitatively compare the performance of the two estimation approaches.
In this section, we formally define the problem set-up in the dynamic pricing setting under consideration. Consider a selling horizon of $T$ time periods, where for each $t = 1, 2, \ldots, T$, the seller offers an itinerary $i$ at price $P_{it}$. Assuming that customers arrive according to a Poisson process and considering an exponential form for the demand response model, the demand and observed prices are modeled as
Here we denote the number of observed purchases for this itinerary in time period $t$ by $Y_{it}$. The vector of controls, denoted by $\boldsymbol{X_{it}}\in \mathbb{R}^m$, includes exogenous features like data related to the specific itinerary derived features or pre-specified transformations of this data (e.g., time seasonality basis functions, constant term etc.), and may also contain lags of $Y_{it}$ and $P_{it}$. $\boldsymbol{W_{it}} \in \mathbb{R}^d$ is the feature set related to the heterogeneity associated with the price-elasticity. $\boldsymbol{W_{it}}$ is based on a sub-vector of the observed data vector $\boldsymbol{X_{it}}$; it includes an intercept term and any pre-specified transformation of raw features. $\boldsymbol{\theta} \in \mathbb{R}^d$ is a vector of parameters associated with price-sensitivity. $\phi(\cdot)$ is a function associated with the total demand volume. To illustrate its role, note that if we set $P_{it} = 0$ in (ref), the purchases will occur according to a Poisson process with mean $\mathtt{exp}(\phi(\boldsymbol{X_{it}})))$. As discussed earlier, the overall demand for a product could be influenced by a variety of factors, therefore $\phi(\cdot)$ could be non-linear, complex and high-dimensional. Given this model, our main objective is to estimate the price-sensitivity parameters $\boldsymbol{\theta}$.
In (ref), we have specified the dependence of price on the set of features $\boldsymbol{X_{it}}$. Typically, with the RM and Pricing systems employed by the firm, the prices are based on marginal opportunity costs (bid prices) computed by solving a finite horizon optimization problem considering demand forecasts, capacities etc., and further refinement based on pricing rules employed by the firm. Therefore, the function $g(\cdot)$ in most practical scenarios will be complex and high-dimensional. Although (ref) is not the main equation of interest, it is useful in elucidating the causal dependence between variables in the model and will also help us later in developing a robust two-stage estimation approach. While our main concern is to estimate the price sensitivity parameters $\boldsymbol{\theta}$, we still need to estimate the additional unknown parameters of the model, $((\phi(\cdot), g(\cdot))$, known in the econometrics parlance as the nuisance parameters. As compared to a fully parametric model such as the Poisson generalized linear model, which itself may lead to biased estimates of $\boldsymbol{\theta}$ if the nuisance (volume) parameter is mis-specified due to linearity restriction, the semi-parametric model in (ref) with a more flexible structure allows one to mitigate this problem while still being interpretable and amenable to efficient estimation approaches for the parameters of interest Robinson1988.
Estimating the parameters of the Poisson semi-parametric model described in (ref) can in practice be challenging. Ideally, one wants $\boldsymbol{X_{it}}$ to incorporate all the possible features and confounders (to control for them), including possible interactions and transformations. To alleviate the challenges related to complex feature engineering, it is beneficial to utilize sophisticated non-parametric ML models in constructing the estimators for the nuisance parameters, $((\phi(\cdot), g(\cdot))$. We next propose two approaches based on these ideas for estimating the parameters of interest.
In the first approach, we model (ref) using a deep neural network and directly estimate the parameters, $\boldsymbol{\theta}$ and $\phi(\cdot)$, by training it on the historical sales transaction data. We would like to leverage the advantages of complex feature representations learned via the deep neural network for estimating the nuisance parameter $\phi(\cdot)$. However, since the dependency of price sensitivity on the features is log-linear, we need to consider a specific neural network architecture for this purpose. The Wide & Deep neural network architecture proposed by Cheng2016 provides a way to accomplish this by jointly training a wide linear part and a deep neural network. Although this architecture was originally proposed for recommender systems, we adapt it for the price-sensitivity estimation problem. Figure (ref) shows this architecture for the Poisson semi-parametric model in (ref).
In this architecture, the wide component is used to model the linear part of the form $P_{it} \boldsymbol{\theta}^T\boldsymbol{W_{it}}$ and the deep part is used to model the nuisance parameter $\phi(\boldsymbol{X_{it}})$. These parts are combined to construct the log rate parameter and fed to a Probabilistic Poisson layer with Poisson negative log likelihood loss for joint training. The estimated parameters of interest, $\hat{\boldsymbol{\theta}}$, can be obtained from the weights of the layer constructing the log-rate parameter.
This method provides a direct way to construct an estimator for the Poisson semi-parametric model and can be implemented relatively easily using packages such as {\tt tensorflow} or {\tt PyTorch}. However, in a practical implementation, this approach can be challenging to work with since the output of the deep part is the nuisance parameter, $\phi(\cdot)$, which is not directly observable. Therefore, we cannot validate how well the deep part is learning this parameter. Moreover, to control for overfitting, a regularization approach (e.g., L1/L2 norm penalization of model parameters, early stopping, drop-out or batch normalization) is often necessary while training the neural network. But regularization may also introduce bias in the estimate of the nuisance parameter. In the absence of properties that make the estimation of price sensitivity parameters robust to errors in the nuisance parameter, any bias in learning the nuisance parameter will lead to bias in the estimate of price-sensitivity parameters $\boldsymbol{\theta}$.
An alternative approach to the Wide & Deep architecture proposed above is to directly use a sophisticated ML estimator (e.g., deep neural networks or gradient boosted trees) for constructing the estimator for $\phi(\cdot)$ in (ref) and estimate the price sensitivity parameters by an iterative scheme that alternates between estimating the parametric and non-parametric parts. We can start by initializing the estimator of $\phi(\cdot)$ to $\hat{\phi}_0(\cdot)$ e.g., by randomly initializing the weights of a neural network. Given the initial estimate of $\hat{\phi}_0(\cdot)$, we can estimate the price-sensitivity parameters, $\hat{\boldsymbol{\theta}}_0$, by maximum likelihood estimation. Given the estimated $\hat{\boldsymbol{\theta}}_0$ we can then create an updated estimate of $\phi(\cdot)$, $\hat{\phi}_1(\cdot)$. This process continues until some convergence criteria is met for the price-sensitivity parameter. However, this approach not only suffers from the regularization related issue described above but due to its iterative nature can also become prohibitively computationally expensive for practical implementation.
To address deficiencies of the estimation approaches discussed so far, we present a two-stage approach in this section. Given the data $\mathcal{D} = (Y, \boldsymbol{X}, P)$, the log-likelihood function for our original model in (ref) can be written as
Let $\boldsymbol{\Theta} \subset \mathbb{R}^d$ be the convex feasible space for price-sensitivity parameters and $\Phi$ be the set of functions of $\boldsymbol{X}$ with finite mean squared error. The true values of the parameters of (ref) are,
where the expectation is taken w.r.t. the probability law governing the data $\mathcal{D} = (Y, \boldsymbol{X}, P)$. Since the log-likelihood function is concave, it satisfies the first order optimality condition given the true parameters $(\boldsymbol{\theta}_0, \phi_0)$
The score function, $\psi(\mathcal{D}; \boldsymbol{\theta}, \phi) = \frac{\partial}{\partial {\boldsymbol{\theta}}}\mathbb{E}[\ell(\mathcal{D}; \boldsymbol{\theta}, \phi)]$, for estimating $\boldsymbol{\theta}_0$ satisfies the Neyman orthogonality condition if the Gateaux derivative of the score function with respect to the the nuisance parameter $\phi$ at $\phi_0$ is 0 chernozhukov2018double, neyman1979c i.e.,
The condition in (ref) intuitively means that the score function for estimating the parameter of interest, $\boldsymbol{\theta}$, is locally insensitive to small perturbations in the nuisance parameter around $\phi_0$. However, the score function, $\psi(\mathcal{D}; \boldsymbol{\theta}, \phi)$, corresponding to our original Poisson semi-parametric model in (ref) doesn't satisfy this orthogonality condition. Thus, we next describe a two-stage estimation method which leads to a reduced form model with score function that is indeed Neyman orthogonal.
In chernozhukov2018double the authors proposed a two-stage estimation approach where in the first stage the nuisance parameters are estimated using sophisticated ML estimators and in the second stage the parameters of interest e.g., the price-sensitivity parameters in our case, are estimated. The second stage score function is constructed so that it satisfies the Neyman orthogonality condition (ref). Most of the works on Neyman orthogonal scores e.g., chernozhukov2018double, semenova2017estimation, consider linear or log-linear models with additive noise. For the Poisson model in our case with non-linear demand response, we construct the Neyman orthogonal scores approximately via the Concentrating-out approach proposed in chernozhukov2018double. This approach works on the principle that maximizing the expected maximum likelihood function with respect to the nuisance parameter, $\phi$, plugging that maximizer back in, and differentiating it with respect to the parameter of interest $\mathbf{\theta}$ produces an orthogonal moment condition.
For this let
be the “concentrated-out” non-parametric part of the model. Note that there is a $\phi_{\boldsymbol{\theta}}$ for any $\boldsymbol{\theta} \in \boldsymbol{\Theta}$.
For our model, this can be derived based on first order conditions as
where the approximation is based on a first order Taylor series expansion of the exponential function around $\mathbb{E}[P|\boldsymbol{X}]$ (details in Section (ref) of the Online Appendix). For the price-sensitivity estimation task we are interested in, this approximation works reasonably well in practice where the term $(\boldsymbol{\theta}^T W\sigma_{\epsilon_{it}})^2$ is small.
Based on the approximate “concentrated-out” non-parametric part, we can construct an alternate score function for estimating $\boldsymbol{\theta}$ which is Neyman orthogonal chernozhukov2018double as
For the Poisson semi-parameteric model, based on (ref) and (ref), the Neyman orthogonal score is given by
Note that the Neyman orthogonal score function for estimating $\boldsymbol{\theta}$ in (ref) corresponds to the reduced form model
Therefore, instead of the original model, we can estimate the price-sensitivity parameters based on the reduced form equation (ref). The subtraction of $\mathbb{E}[P_{it}|\boldsymbol{X_{it}}]$ from $P_{it}$ in (ref) changes the price-related variable to have a “local" interpretation, in which the parametric part locally governs changes in demand w.r.t. changes in price around the expected value for price given $\boldsymbol{X_{it}}$.
Working with the reduced form version of the demand response model is advantageous for several reasons. Firstly, instead of the nuisance parameters $\phi(\cdot)$ and $g(\cdot)$ in (ref), the nuisance parameters in (ref) are expected purchases $\mathbb{E}[{Y}_{it}|\boldsymbol{X}_{it}]$ and expected price $\mathbb{E}[{P}_{it}|\boldsymbol{X}_{it}]$. These are directly observable quantities so we can deploy sophisticated ML techniques for constructing estimators $\hat{Y}_{it}$ and $\hat{P}_{it}$ for $\mathbb{E}[{Y}_{it}|\boldsymbol{X}_{it}]$ and $\mathbb{E}[{P}_{it}|\boldsymbol{X}_{it}]$; validation of their goodness of fit becomes much easier, enabling us to use standard techniques for model selection and hyperparameter tuning. Furthermore, learning $\hat{P}_{it}$ is akin to learning the firm's historical pricing policy which can be useful in constructing reasonable price bounds for an automated dynamic pricing system in a practical deployment. Secondly, given estimates $\hat{P}_{it}$ and $\hat{Y}_{it}$, the model in (ref) is now a fully parametric Poisson generalized linear model, therefore it is easy to use batch or online estimators for estimating the price sensitivity parameters $\boldsymbol{\theta}$. Finally, and most importantly, by approximately orthogonalizing the price with respect to $\boldsymbol{X}$ and removing the effect of confounding by subtracting $\hat{P}_{it}$ from $P_{it}$ in (ref), the reduced form model provides robustness in the subsequent estimation of parameters $\boldsymbol{\theta}$ against biases in the estimators for nuisance parameters chernozhukov2018double.
Our overall scheme for constructing the estimator of price-sensitivity parameters $\boldsymbol{\theta}$ based on the reduced form (ref) comprises of two stages. In the first-stage, we create two estimators, $\hat{P}_{it}$ and $\hat{Y}_{it}$, for $g(\boldsymbol{X_{it}}) = \mathbb{E}[P_{it}|\boldsymbol{X_{it}}]$ and $\mathbb{E}[Y_{it}|\boldsymbol{X_{it}}]$, respectively, using machine learning methods. For $\hat{P}_{it}$ the response variable is the observed price and for $\hat{Y}_{it}$ the response variable is the observed bookings. Given the estimates from the first stage, in the second stage we estimate parameters of a Poisson generalized linear model. The use of non-parametric ML-based methods for the estimators of $\hat{P}_{it}$ and $\hat{Y}_{it}$ allows us to bring more complex representation of the feature space $\boldsymbol{X_{it}}$ into the model while still being able to perform the second-stage estimation via a tractable and interpretable Poisson generalized linear model form.
In this work, we use neural networks for constructing the first-stage estimators and Bayesian dynamic generalized linear models (DGLM) for the second stage estimation as will be discussed in the following sections. The use of Bayesian methods in the second stage is mainly done to get online estimators of $\boldsymbol{\theta}$ which are convenient for practical applications. Any concern regarding the regularization bias due to implicit regularization in Bayesian methods can be mitigated by considering prior distributions with high variance. Furthermore, the Bayesian approach also provides a principled way to quantify uncertainty on the estimate of price-sensitivity parameters which is not only useful for validating the model but can also drive reinforcement learning-based approaches for price optimization. Similar to chernozhukov2018double, semenova2017estimation, athey2019generalized, and mackey2018orthogonal, we use cross-fitting approach in our two-stage scheme to prevent overfitting, where the data is split into $K$-folds and for each fold the first stage predictions used in the second stage estimation are from first-stage models trained on other folds. As shown in chernozhukov2018double prevention of over-fitting in the construction of first stage estimators for $\mathbb{E}[Y_{it}|\boldsymbol{X_{it}}]$ and $\mathbb{E}[P_{it}|\boldsymbol{X_{it}}]$ plays a key role in getting $\sqrt{n}$-consistent two-stage estimators for $\boldsymbol{\theta}$.
In (ref), we assume a parametric form for treatment heterogeneity because of its interpretability. In practice, a domain expert may choose $\boldsymbol{W_{it}}$ considering the trade-off between flexibility and estimation precision. $\boldsymbol{W_{it}}$ can, for example, be a vector of binary variables encoding categorical features and other continuous features such as basis functions to capture seasonality etc. In our simulation experiments, we empirically see that our method provides elasticity estimates that are robust to minor misspecification in $\boldsymbol{W_{it}}$ when compared to simpler approaches. We will refer to this two-stage estimation approach as Two-stage-CO in the rest of this article.
For real-world applications, we use neural networks to construct the estimator for $\mathbb{E}[P_{it}|\boldsymbol{X_{it}}]$ in the first stage with the mean-squared error loss. Although in principle we can use other function approximators such as gradient boosting, neural networks can provide advantages beyond just the prediction task. One such advantage is the ability to encode categorical variables in our feature set via learnable embedding vectors. We can cluster categorical variables based on the closeness (euclidean norm) of embedding vectors e.g., derive groups or segments related to Origin-Destination pairs, Days of Weeks etc., with similar pricing behavior. Such inferred closeness of categories contains information about expert knowledge encoded in the pricing system, which then can be used to inform the choice of the structure of $\boldsymbol{W}$ and/or segmentation to help with data sparsity.
When learning the first stage model for price, we leverage all the available data from different markets (Origin-Destination pairs) and add categorical variables representing each market to the feature set. This also requires careful scaling of the prices in different markets to a similar range so that they have a balanced influence on the training. By pooling all the data, we increase the sample size and also effectively add an implicit regularization to learn a more accurate model compared with training on each market separately. In addition, we embed the categorical features in a low-dimensional continuous space, where the embedding is learned jointly with other parameters of the neural network.
We also estimate $\mathbb{E}[Y_{it}|\boldsymbol{X_{it}}]$ using neural networks in the first stage. However, as the generating process for purchases is assumed be a count process, we train the model with a Poisson log-likelihood loss. Note that the choice of this loss function is still an approximation as although $Y_{it}|\boldsymbol{X_{it}}, P_{it}$ has a Poisson distribution, the distribution of $Y_{it}|\boldsymbol{X_{it}}$ is not Poisson. This approximation is mainly driven by the ease of implementing neural networks using standard packages such as Pytorch and tensorflow. Our numerical results in Sections (ref) and (ref) show that despite this approximation, the two-stage estimation approach performs well.
Given the predictions of the first-stage models, in the second stage, we adopt a Bayesian dynamic generalized linear model (DGLM) west2006bayesian,berry2020bayesian for estimating the price-sensitivity parameter $\mathbf{\theta}$ in (ref). DGLMs are based on the Bayesian state-space models with a sequential and efficient update procedure, which makes them suitable for applications with streaming data. Additionally, the Bayesian framework provides a principled way for quantifying uncertainty of parameter estimates, which is helpful for making more informed pricing decisions and utilizing RL techniques. Based on the estimates from the first stage, the Bayesian model for the second stage for (ref) is:
Here, $\hat{Y}_{it}$ is the cross-fitted prediction of purchases for the specific itinerary $i$ from the first stage model and $\pi(\boldsymbol{\theta})$ is the prior distribution over the parameters. The inference in the second stage is done independently for each market on all itineraries belonging to that market. The inference procedure of the second stage can be completely parallelized for efficiency, which is extremely relevant for an implementation in a real-life system. For each market, $\boldsymbol{W}_{it}$ is concatenation of an intercept, continuous features (e.g., Fourier basis for departure date seasonality, Splines related to days before departure, departure time etc) and one-hot encoded categorical variables (e.g., Day of Week, Point of Sale, Number of Connections etc.) as shown in (ref).
Thus, by modeling heterogeneous price-sensitivity through an additive structure, the parameter $\boldsymbol{\theta}$ in (ref) is estimated for each market (origin-destination pair) using data for all its departure dates present in the historical transactions.
The update procedure in DGLMs is sequential and can be performed immediately after each new observation (daily bookings and price recorded for a market). The overall procedure is based on the combination of linear Bayes' theory goldstein2007bayes and a variational Bayes approach. We describe the inference procedure in detail in Section (ref) of the Online Appendix, where we also introduce a Laplace approximation based approach to further speed up the procedures of berry2019bayesian and west2006bayesian. Through this inference procedure, we are able to update the posterior means and covariances of the parameters without requiring any assumptions on their exact distribution.
A desirable property of the dynamic regression model is that it can be incrementally updated in an efficient manner in between the complete retraining of the first and second stage models. In our application, the updating schedule of the first stage and second stage models are different. Since the airlines pricing policies and the underlying model of latent variables are expected to be more stable, we update the first stage models, which have higher computational cost, with lower frequency compared with the second stage model. The second stage model is incrementally updated after each observation to have the most updated belief on $\boldsymbol{\theta}$.
An alternative approach to construct Neyman orthogonal scores for semi-parametric generalized linear models has recently been proposed in nekipelov2022regularised. In this approach, the authors first partial out Robinson1988 the effect of controls ($\boldsymbol{X}$) from treatment ($P$) inside the link function's argument of the GLM. Specifically, for the case of Poisson semi-parametric model (ref), in this approach one starts with the following partialled-out form for the original model (ref):
where $q(\boldsymbol{X}_{it}) := \mathbb{E}[P_{it}|\boldsymbol{X}_{it}]\boldsymbol{\theta}_0^T\boldsymbol{W_{it}} + \phi_0(\boldsymbol{X}_{it}) = \mathbb{E}[\mathtt{log}(\mathbb{E}[Y_{it}|\boldsymbol{X_{it}},P_{it}])|\boldsymbol{X}_{it}]$. The Neyman orthogonal score for the Poisson semi-parametric model based on nekipelov2022regularised is:
which corresponds to a weighted maximum likelihood based score function for (ref) with the weights being $\frac{1}{\mathbb{E}[Y|\boldsymbol{X},P]}$.
To utilize this score function, we first need to construct (ML) estimators for $\mathbb{E}[P_{it}|\boldsymbol{X}_{it}]$ and $\mathbb{E}[Y_{it}|\boldsymbol{X}_{it},P_{it}]$. The predictions from the initial estimator for $\mathbb{E}[Y_{it}|\boldsymbol{X}_{it},P_{it}]$, are used as weights and also as response variables to form an estimator for $q(\boldsymbol{X}_{it})$, i.e. an estimator for $\mathbb{E}[\mathtt{log}(\mathbb{E}[Y_{it}|\boldsymbol{X_{it}},P_{it}])|\boldsymbol{X}_{it}]$. After constructing these estimators, the parameter of interest, $\boldsymbol{\theta}$, is found based on the score function in (ref). Given the first-stage estimators, the solution to (ref) can be easily found by fitting a regular weighted Poisson generalized linear model. We will refer to this partialled-out approach as Two-stage-PO and the concentrated-out approach in Section (ref) as Two-stage-CO. It is worth noting that the partialling-out approach involves building (ML) estimators for $\mathbb{E}[P_{it}|\boldsymbol{X}_{it}]$, $\mathbb{E}[Y_{it}|\boldsymbol{X}_{it},P_{it}]$, and $q(\boldsymbol{X}_{it})$, while the Two-stage-CO approach only needs (ML) estimators for $\mathbb{E}[P_{it}|\boldsymbol{X}_{it}]$ and $\mathbb{E}[Y_{it}|\boldsymbol{X}_{it}]$. Thus, the Two-stage-CO method can be significantly better in terms of computation time. Moreover, the inclusion of an additional estimator adds to the complexity of the process and requires increased effort in validating the goodness of fit of these estimators. We will compare the estimation performance of Two-stage-CO and Two-stage-PO via Simulation studies in Section (ref).
In this section we describe how to compute dynamic pricing policy given the estimates of the price-sensitivity parameters. The objective of the firm is to maximize the expected margin contribution given the associated per unit cost $c$ from selling the product to a request during a shopping session. As noted earlier, for the airline pricing problem, the cost is the marginal opportunity cost or bid price associated with the seat inventory requested in the itinerary. Here we assume that the applicable marginal opportunity cost has already been computed by the airline's RM system by solving a separate finite horizon optimization problem and is available at the time of the request.
Let $\boldsymbol{\mu}_{\theta}$ and $\Sigma_{\theta}$ be the current estimate of mean and covariance associated with the price-sensitivity parameters $\boldsymbol{\theta}$ at the time of request, and $\mathcal{P}_{it} := p^{lb}_{it}\leq p \leq p^{ub}_{it}$ be the feasible space of prices for the request. Let $R_{it}(p; c, \boldsymbol{X}_{it}, \boldsymbol{\mu}_{\theta},\Sigma_{\theta})$ denote the expected margin contribution from an itinerary/product $i$ requested at time $t$ with associated feature vector $\boldsymbol{X}_{it}$ and cost $c$. The optimal price for the request is given by
In Section (ref) we discuss how we compute the dynamic prices in our experiments.
Given the posterior distribution of price sensitivity parameter, $\pi^{*}(\boldsymbol{\theta})$, the expected margin contribution for the Poisson demand response model in (ref) is:
Finding the price based on maximizing the above expected margin contribution is akin to using a greedy policy for pricing. Since the Bayesian DGLM-based update procedure in (ref) does not make distributional assumptions on $\pi^{*}(\boldsymbol{\theta})$ and only yields the posterior mean and covariance of $\boldsymbol{\theta}$, a closed-form expression of (ref) is not available, and one can either resort to a plug-in approach using point estimates of the parameters or an approximation of the expectation. We discuss the plug-in approach that we use in our experiments in the following and two more approaches to approximate (ref) in Section (ref) of the Online Appendix.
Using the posterior mean of the parameters, which is the optimal mean-squared error estimator, we can form a plug-in estimator of (ref), where the corresponding margin contribution function is
The utility function in (ref) results in a convex optimization problem that in fact has a closed-form solution
Note that if we employ frequentist methods for estimating $\boldsymbol{\theta}$ e.g., the direct estimation approach in Section (ref), we can use the same formula in (ref) to compute the optimal price which maximizes the margin contribution by replacing $\boldsymbol{\mu}_{\theta}$ with the frequentist estimate $\hat{\boldsymbol{\theta}}$.
Some other methods for constructing pricing policies such as Bayes greedy and Reinforcement Learning are discussed further in Section (ref) of the Online Appendix.
In this section we present results from computational experiments on randomly generated data where we know the true ground demand model, i.e., we know the underlying price-demand relationship. The goal of this study is to analyze and compare the accuracy of the proposed estimators of price-sensitivity parameters based on simulated data.
We first perform a comparison on a simple setup with randomly generated data where we set the ground truth of the parameters. By comparing the estimates with the true parameters, we can evaluate the performance of our proposed estimation methods. For this simulation, the exogenous features, $\boldsymbol{X} \in \mathbb{R}^{10}$, are normally distributed, $\boldsymbol{X}_i \sim \mathcal{N}(0,\Sigma(\rho))$, where $\Sigma_{jk}(\rho)=\rho^{|j-k|}$ and $\rho$ is set to 0.5. The dependence of price on the features is assumed to be linear and the demand data are generated according to a Poisson model:
Here, $\tau=1.2$, $\boldsymbol{\eta} \in \mathbb{R}^{14}$ with all elements equal to 0.1, and $\boldsymbol{\gamma} \in \mathbb{R}^{10}$ with all elements equal to 3. As can be seen, the nuisance $\phi(\boldsymbol{X_i})$ is a simple polynomial of degree 2. The vector of parameters of interest, $\boldsymbol{\theta}^T=[\theta_0,\theta_1,\theta_2,\theta_3,\theta_4]$, is set as [-0.02,-0.005,-0.005,-0.005,-0.005].
We compare the performance of the direct estimation approach based on Wide and Deep architecture presented in Section (ref), and the two-stage estimation approaches discussed in Section (ref). For this example, the deep part of the Wide and Deep consists of one hidden layer with 50 units and ReLU activation followed by a layer with a single unit representing the nuisance parameter. Additionally, we use L2 regularization on hidden layer weights of the deep part, and in each run, we do early stopping and restore the best performing model state on a validation set, formed by randomly setting aside 20% of training data. These are all techniques commonly used when training neural networks on real data which can introduce some level of bias in parameter estimates. The model is trained using the Adam optimizer kingma2015adam.
For the proposed two-stage methods, the reduced form model corresponding to the Two-stage-CO approach in Section (ref) is:
We use ridge regression for estimating $\mathbb{E}[P_{i}|\boldsymbol{X_{i}}]$ in the first stage. The first stage estimator for $\mathbb{E}[Y_{i}|\boldsymbol{X_{i}}]$ is constructed using a random forest regressor with mean-squared-error loss. All first stage model training and predictions are carried out by 5-fold cross-fitting. In this toy example, to isolate the effect of the reduced form model on estimation accuracy, we employ a frequentist Maximum Likelihood Estimation (MLE) approach for the Poisson GLM in the second stage estimation of the parameter of interest, $\boldsymbol{\theta}$.
For the alternative Two-stage-PO estimation approach discussed in Section (ref), in the first stage, we first fit $\mathbb{E}[Y_{i}|\boldsymbol{X_{i}},P_i]$ with a random forest regressor with Poisson deviance loss, and then construct an estimator for $\mathbb{E}[\mathtt{log}(\mathbb{E}[Y_{i}|\boldsymbol{X_{i}},P_i])|\boldsymbol{X_{i}}]$ by training a random forest regressor with mean-squared-error loss on $\mathtt{log}$ of predictions of the previous random forest. Note that these steps are performed with 5-fold cross-fitting. To estimate the parameter of interest, in the second stage we use a frequentist weighted MLE approach for the Poisson GLM where samples are given weights inversely proportional to their first stage prediction of $\mathbb{E}[Y_{i}|\boldsymbol{X_{i}},P_i]$.
We generate 10,000 samples from the assumed data generation model in each simulation run, and test the estimation accuracy of different methods on the generated data. The procedure of data generation and model testing is repeated ten times and the average results are shown in Table (ref).
The two-stage methods result in significant improvement over the direct estimation procedure. Moreover, Two-stage-CO and Two-stage-PO show very close performance and the difference between them is not statistically significant.
To generate the sales transactions mimicking the complexities involved in real-life airline pricing, we model the interactions between RM and pricing with finite seat capacity. We set-up the simulation experiment considering a single-leg with a seat capacity of 100. The booking horizon for each departure date consists of 365 days grouped into 10 contiguous time periods called time frames (TF), such that the demand arrival rate within each TF is constant but varies across TFs. TF 0 is the furthest away from departure of a flight and TF 9 is closest to departure. The demand arrival is Poisson distributed with a rate that varies by week of year (WOY), booking period (TF), days of week (DOW) and point of sale (POS). Moreover, customer willingness-to-pay is exponentially distributed ($\mathtt{WTP}|\theta \sim \mathtt{exp(-\frac{1}{\theta})}$) and the price-sensitivity parameters for customers ($\theta$) vary by POS and TF. With this setup, the ground demand response model governing demand on day $t$ of the booking horizon for departure date $i$ is
The true price-sensitivity parameters are shown in Table (ref). Note that as discussed previously, for the exponential demand response function, the optimal price given price sensitivity parameter $\theta$ and cost $c$ is given by
For a given POS, TF combination, $\alpha = -\frac{1}{\theta}$ denotes the mean willingness-to-pay, which is also the optimal price when the cost is zero. These provide a more intuitive explanation of the price sensitivity parameters and therefore we have included them in the Table (ref).
The RM system computes the marginal opportunity cost of seats (bid prices) considering the known demand arrival rates and price-sensitivity parameters using dynamic programming and the pricing system computes the final prices based on the relevant bid price and true price-sensitivity parameters i.e., the optimal pricing policy. To mimic idiosyncratic noise in prices observed in real life, we add normally distributed noise to final prices with mean zero and standard deviation of 20. These prices are then offered to a stream of arriving customers and the customer accepts or rejects the price based on their sampled willingness-to-pay ($\mathtt{WTP}|\theta \sim \mathtt{exp(-\frac{1}{\theta})}$). The number of customers arriving in any time unit and their WTP are sampled according to the encoded ground truth price-demand model. This process creates the historical sales transaction data which includes daily bookings, daily average price and additional features useful for estimation namely, $\mathtt{WOY, DOW, POS, TF}$. Given this scheme, the price generation process in the simulation is:
where $v_{it}(x)$ denotes the optimal opportunity cost (bid price) computed via dynamic programming based optimization for the departure date $i$ on day $t$ given that we have $s_{it}$ seats remaining. Given the complex dependence of price on bid price, reservation process and price elasticity, the data generation process for bookings and prices resembles the model described in (ref) and (ref).
Using this simulation setup, we generate data corresponding to two years worth of daily flight departures (730 individual flights). Further analysis on the data generated under this simulation set-up can be found in the Online Appendix, Section (ref).
For the purpose of estimating the price-sensitivity parameters, we consider the following model for the log-rate of the Poisson demand response model
Note that the ground demand response model used for the data generating process has a weak interaction in the price sensitivity parameters between POS and TF. Given some of the sparser entities, estimation benefits from a simpler linearly additive specification of price effects across POS and TF as shown in (ref). Moreover, minor misspecifications of this type are common in practice, so the model in (ref) will allow us to validate the robustness of the estimation method better. We use the direct neural network based technique and the two-stage approach to estimate the parameters of interest, $\boldsymbol{\theta}$, for the above model next.
We first perform the direct estimation method proposed in Section (ref). The features used in the model, $\boldsymbol{X}_{it}$, include one-hot encoded categorical variables for POS, TF, and DOW, and four continuous features representing the yearly seasonality, which are Fourier series based on WOY with the first two frequencies and 52 weeks period. The deep part consists of a hidden layer with 50 units and ReLU activation followed by a layer with a single unit. The layer constructing the Poisson log rate forms $P_{it}(\sum_{pos=1}^{\#POS-1}\boldsymbol{\theta}_{pos}\mathbbm{1}_{\{POS=pos\}}+ \sum_{tf=0}^{\#TF}\boldsymbol{\theta}_{tf}\mathbbm{1}_{\{TF=tf\}}) + \phi(\boldsymbol{X}_{it})$. In each run, we randomly set aside 15% of flights as validation set and train the model using the Adam optimizer kingma2015adam with early stopping.
As detailed in Section (ref), we construct two-stage estimators for the price sensitivity parameters $\boldsymbol{\theta}$ using both Two-stage-CO and Two-stage-PO methods.
For the Two-stage-CO approach detailed in Section (ref), we use the following reduced form model based on (ref)
The first-stage estimators for $\mathbb{E}[{P}_{it}|\boldsymbol{X}_{it}]$ and $\mathbb{E}[{Y}_{it}|\boldsymbol{X}_{it}]$ are constructed using a random forest regressor where the number of decision trees to be ensembled are fixed to 100 and use a mean squared error loss function. The features are the same as the ones used in the direct method explained in the previous Section. Moreover, to avoid overfitting, we use a cross-fitting approach and split the data into 5-folds with random shuffling. When estimating the second-stage model on each fold, we use the predictions from the first-stage models trained on other folds.
For the second-stage estimator we use the Bayesian Dynamic Generalized Linear Model approach described in Section (ref). For the model parameters, $\boldsymbol{\theta}$, we set the prior mean vector, $\boldsymbol{\mu}_{\boldsymbol{\theta}} = \boldsymbol{0}$, and a diagonal covariance matrix with all variance values set to 10 and the discount factors $\delta^{\boldsymbol{\theta}} =1.0$.
The alternative Two-stage-PO approach detailed in Section (ref) we require three first-stage estimators. For the estimator of $\mathbb{E}[Y_{i}|\boldsymbol{X_{i}},P_i]$ we utilize a random forest regressor with Poisson deviance loss, and then construct an estimator for $\mathbb{E}[\mathtt{log}(\mathbb{E}[Y_{i}|\boldsymbol{X_{i}},P_i])|\boldsymbol{X_{i}}]$ by training another random forest regressor with mean-squared-error loss on $\mathtt{log}$ of predictions of the previous random forest. In addition the estimator for $\mathbb{E}[{P}_{it}|\boldsymbol{X}_{it}]$ is constructed in a similar manner to that in the Two-stage-PO approach. In all the random forest regressors, the number of decision trees to be ensembled are fixed to 100. The second-stage estimation for the price-sensitivity parameters in the Two-stage-PO approach is performed using a frequentist maximum-likelihood estimation. We use a frequentist approach here since for the score function in (ref) we need to perform a weighted maximum likelihood estimation using the estimates of $\frac{1}{\mathbb{E}[Y_{i}|\boldsymbol{X_{i}},P_i]}$ as weights. Therefore, this weighting scheme cannot be implemented in a meaningful way via Bayesian modeling.
The final values of estimated price sensitivity parameters $\hat{\alpha} = -\frac{1}{\hat{\theta}}$ after updating over 104 week's worth of data for both direct and two-stage methods are shown in Table (ref). The mean absolute percentage errors (MAPE) between the estimated and true price sensitivity parameters are also shown in Table (ref). We also compute the weighted MAPE (wMAPE) metric where the weights are normalized to 1 based on relative bookings received for various POS/TF combinations i.e., we give higher weight to POS/TF combinations with more bookings. Since the direct approach's performance may change based on the train-validation split and different random initialization of the model, we perform 10 runs and for each run calculate the MAPE and wMAPE and report the means in Table (ref). Note that the individual numbers for the direct approach estimates in Table (ref) are from one run that had performance close to the mean (MAPE=23.814% and wMAPE=0.786%) and are used as a representative for illustration purposes.
We see that both the direct approach via wide and deep neural network and the two-stage approaches are able to estimate price-sensitivity parameters in a reasonable range and capture general trends like willingness-to-pay increasing as we get closer to departure and that for a given TF, WTP for POS 1 are higher than that for POS 0. However, the two-stage approaches have lower MAPE errors of $4.16\%$ and $17.36\%$ for Two-stage-CO and Two-stage-PO respectively, as compared to the Wide and Deep approach with an average MAPE error of $24.876 \%$. From the wMAPE metric perspective, the Two-stage-CO approach is substantially better than the Two-stage-PO and the Direct approach which have very similar wMAPE error. These results also show the robustness of the two-stage approach to the minor mis-specification in the linear treatment effect model.
We note that the Two-stage-CO method provides substantial improvement in the estimates of price sensitivity parameters as compared to other approaches. The estimation for the Two-stage-CO was performed using Bayesian method while that for Two-stage-PO was done via frequentist approach. To rule out the possibility of difference in second stage estimation methodology leading to difference in performance of these two methods, we also performed the estimation for the Two-stage-CO via the same frequentist Poisson GLM approach. We can see from Table (ref) provided in Section (ref) of the Online Appendix that the performance of Bayesian and frequentist estimates for the Two-stage-CO is very close. Therefore, the difference in the performance of these methods cannot be attributed just to the Bayesian approach utilized for the second stage.
Although the Two-stage-CO and Two-stage-PO approaches performed similarly in simpler simulation setting, in more complex settings the added estimation required in the first-stage for the Two-stage-PO may compound errors and lead to deterioration of estimation performance. These results indicate that the Two-stage-CO is not only computationally efficient but also leads to excellent estimation performance in both simple and complex simulation settings. Therefore, for the numerical study in the next section we choose to perform analysis only on the Direct and the Two-stage-CO approach.
The Bayesian DGLM model in the Two-stage-CO approach updates the parameter estimates in an online fashion, Figure (ref) shows the evolution of the posterior means of price sensitivity parameters ($\hat{\alpha} = \frac{1}{\hat{\mu}_{\theta}}$) for various (POS, TF) combinations. We see that most of the parameters converge to a value close to true parameter value within 20 weeks worth of data. The sparser (POS, TF) combinations e.g., TF 8 take a longer time to converge.
In this section, we present results of the proposed methodologies on an anonymized data set comprising of airline sales transactions using an offline estimation study. For real data, since the true price-sensitivity parameters are unknown, it is difficult to validate the quality of the estimators in a manner similar to the simulation study. Instead, we perform a qualitative study to determine the reasonability of the recommended optimal price based on the estimated price-sensitivity parameters by comparing these to the historical prices in the sales transaction data. We also discuss reasonability of trends and patterns observed in the price-sensitivity parameters vis-a-vis expected trends.
The data comprises of historical transactions for a period of 100 weeks for four markets (Origin-Destination pair). The transaction data is for economy compartment and includes daily bookings and daily average price for each departure date for a 180 day booking horizon.
In building the model for price sensitivity estimation, we considered several features including origin and destination, departure time, travel time. Furthermore, departure date related features like yearly seasonality via Fourier series representations based on departure day of year, departure day of week, holidays/special events and booking date related features like seasonality via Fourier series representations based on days before departure were also included. For estimating price sensitivity on the real data, we consider price sensitivity to be differentiated by POS, DOW, and Days before Departure (DBD) for each Origin-Destination pair. We use continuous basis functions for representing price sensitivity variation over DBDs.
For direct estimation by the wide and deep method, we use a three layer architecture for the deep part. The features used in the model include the categorical features, and continuous features representing seasonality mentioned above. The wide and deep architecture is trained for each Origin-Destination separately.
We use neural networks with five layers for the first-stage prediction of price ($\hat{P}_{it}$) and demand ($\hat{Y}_{it}$). The models are trained via a cross-fitting approach, where the flights are split into 10-folds. We observed that by pooling data from other origin-destinations to learn a model for the conditional expected price, the prediction accuracy was improved by 5% for the four origin-destinations of interest, compared to using data only from one origin-destination.
For the second-stage estimator we use the Bayesian Dynamic Generalized Linear Model approach described in Section (ref). For the model parameters, $\boldsymbol{\theta}$, we set the prior mean vector, $\boldsymbol{\mu}_{\boldsymbol{\theta}} = \boldsymbol{0}$, and a diagonal covariance matrix with all variance values set to 2 and the discount factors $\delta^{\boldsymbol{\theta}} = 0.9995$.
Unlike the simulation setting, in the real data setting, we cannot compare the estimated price sensitivity parameters to some ground truth. Nevertheless, we can still make some qualitative statements by comparing the price recommended under various estimation schemes and the historical prices. Figure (ref) shows the recommended price under Direct and Two-stage estimation approach in solid lines and historical prices in dashed lines over the booking horizon (days before departure) for four airline markets.
The historical and recommended prices are averaged over other dimensions like departure dates. Note that, if the estimate of the price sensitivity parameter is $\hat{\alpha} = -\frac{1}{\hat{\mu}_{\theta}^T \mathbf{W}}$, then the optimal price given cost $c$ is $c +\hat{\alpha}$. The appropriate cost to use in our setting is the marginal opportunity cost or the bid price computed by the RM system for the inventory position at the time of the booking request. We use the historical marginal opportunity costs stamped in the booking transactions for computing the recommended price corresponding to the price sensitivities estimated with Direct Wide and Deep Neural Network and the Two-stage-CO methods.
The left panel in Figure (ref) shows the evolution of price recommendations based on price-sensitivities estimated using direct estimation approach. We see from these figures, that the recommended prices from the direct approach increase as we go closer to departure on all markets. However, we observe that the recommended prices are much higher than the historical prices for all the markets. Although we do not expect the prices in the historical transactions to be optimal, if a method generates prices considerably higher or lower than historical prices, it is very likely that the price sensitivity estimates are biased.
The right panel in Figure (ref) shows the price recommendations based on price-sensitivities estimated using the two-stage-CO approach. We see that similar to the direct estimation approach, the recommended prices are increasing as we go closer to departure. However, unlike the direct estimation approach, the aggressiveness of this increase is much lower and more in line with historical prices. We also see that the price trend across different days of the week (DOW) is similar to that seen in the historical prices.
Figure (ref) shows the distribution of the difference between the recommended and historical prices (positive values indicate that the recommended price is higher than the historical price) for the Direct and Two-stage-CO approach. The distributions are plotted differently corresponding to transactions with positive bookings and those with zero bookings. The prices recommended by the Direct estimation method tend to be much higher than the historical prices (mode occurs at more positive value) for both booking and non-booking related transactions. The prices recommended by the Two-stage-CO approach tend to be closer to historical prices (mode is closer to zero). As with the Direct approach, for most of the markets, the mode of the empirical distribution for the Two-stage-CO method also occurs at slightly positive values indicating that the willingness-to-pay may be higher than that used during the historical period. Although qualitative in nature, these results do provide some evidence that the over-estimation of price-sensitivity parameters observed in the direct estimation approach has been corrected to a large extent by the two-stage estimation approach.
As discussed earlier, one of the motivations for keeping the specification of price sensitivity parameters interpretable is that we can easily visualize how these parameters vary across various factors such as departure time, weak of year, days to departure etc. Figure (ref) shows this analysis for Market 1. We can observe the following behaviors for this market:
These observations are aligned with typical behavior on this market and generates further confidence in the estimates of price-sensitivity parameters from the Two-stage-CO methodology.
In this paper, we introduced a Poisson semi-parametric demand response model with feature-dependent price-sensitivity to enable automated dynamic pricing. This model allows one to construct a more flexible model than a fully parametric one. In particular, we allow for the nuisance part to use modern ML techniques while keeping a log-linear specification for the effect of price on demand. After first proposing a direct modeling and estimation approach for this model via a deep neural network architecture, we then propose a two-stage estimation approach for price-sensitivity estimation based on ideas related to double/debiased ML approach and Neyman orthogonalization of score functions. This approach is advantageous not only from the perspective of practical implementation but also endows robustness to the estimates of price-sensitivity from biases in the estimation of the nuisance part. The second-stage inference of parameters driven by a variational Bayes approach leads to a computationally efficient sequential learning scheme and also enables the implementation of policies for smart price-experimentation based on RL techniques. The numerical studies performed based on simulated data as well a real data demonstrate the superior performance of the two-stage estimation approach as compared to the direct estimation approach.
Although this work was applied to the setting of Airline pricing, we believe that the model and associated estimation techniques discussed in this paper are relevant for retail pricing problem in many industries. Furthermore, the general philosophy of combining traditional fully parametric models with more sophisticated ML techniques e.g., via semi-parametric models, to construct powerful frameworks for decision problems will be helpful beyond the pricing setting discussed in this work. In contrast, for thosesettings where interpretability is not an important concern, extension of the two-stage estimation techniques to fully non-parametric models for dynamic pricing is a promising direction for future research.