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.
38,927 characters · 17 sections · 35 citation commands
Assessing the effect of advertising expenditures upon sales: a Bayesian structural time series model
\captionsetup{width=0.8\textwidth}
\markboth
{[email removed]}{1}
{[email removed]}{2}
{[email removed]}{3}
{[email removed]}{{1,2,4}}
{Institute of Mathematical Sciences} {ICMAT-CSIC}
{Facultad de Ciencias F\'isicas}{Universidad Complutense de Madrid}
{ETSIN}{Universidad Polit\'ecnica de Madrid} {Department of Computer Science, School of Engineering}{Universidad de C\'adiz}
It is widely acknowledged that a firm's expenditure on advertising has a positive effect on sales assmus1984advertising, tellis2007advertising, luo2012does, wiesel2011practice. However, the exact relationship between them remains a moot point, see e.g. tellis2009generalizations for a broad survey. Since Dorfman and Steiner's dorfman1954optimal seminal work\,---one of the first formal theories of optimal monopoly advertising---\,several models have been proposed to pinpoint this relationship, although consensus on the best approach has not been reached yet. Two diverging model-building schools seem to dominate the marketing literature little1979aggregate: a priori models that rely heavily on intuition and are derived from general principles, although usually with a practical implementation on mind (nerlove1962optimal, or vidale1957operations and little1975brandaid inter alia) and statistical or econometric models, which usually start from a specific dataset to be modelled (see assmus1984advertising for a review). In this work we will mostly rely on the first type of models, viz. that of Nerlove-Arrow nerlove1962optimal, which extends Dorfman and Steiner to a dynamic setting bagwell2007economic and adapts seamlessly to the \emph{state-space} or \emph{structural} time series approach.
Bayesian structural time series models scott2014predicting, in turn, have positioned themselves in the past few years as very effective tools not only for analysing marketing time-series, but also to throw light into more uncertain terrains like causal impacts, incorporating a priori information into the model, accommodating multiple sources of variations or supporting variable selection. Although the origins of this formalism can be traced back to the 1950's in engineering problems of filtering, smoothing and forecasting, first with Wiener wiener1949extrapolation and specially with Kalman kalman1960new, these problems can also be understood from the perspective of estimation in which a vector valued time series $\{ X_0, X_1, X_2, \ldots\}$ that we wish to estimate (the latent or hidden states) is observed through a series of noisy measurements $\{ Y_0, Y_1, Y_2, \ldots\}$. This means that, in the Bayesian sense, we want to compute the joint posterior distribution of all the states given all the measurements sarkka2013bayesian. The ever-growing computing power and release of several programming libraries in the last few years like petris2010r, scott2016bsts have in part alleviated the difficulties in the mathematical underpinning and computer implementation that this formalism suffers, making these methods broadly known and used. This family of models have been used successfully to model financial time series data doi:10.1002/asmb.428, infer causal impact of marketing campaigns brodersen2015inferring, select variables and nowcast consumer sentiment scott2015bayesian, or for predicting other economic time series models like unemployment scott2014predicting.
In this paper we use the formalism of Bayesian structural time-series models to formulate a robust model that links advertising expenditures with weekly sales. Due to the flexibility and modularity of the model, it will be well suited to generalization to various markets or situations. Its Bayesian nature also adapts smoothly to the issue of introducing outside or a priori information, which can be updated to the posterior distribution of the estimated parameters. The formulation of the model allows for non-gaussian innovations of the process, which will take care of the heavy-tailedness of the distribution of sales increments. We also discuss how the forecasts produced by this model can help the manager into allocating the advertising budget. The decision space is reduced to a one dimensional curve of Pareto optimal strategies for the two moments of the forecast distribution: expected return and variance.
The paper is organized as follows: after a brief review of the most usual marketing models and the formalism of structural time series in section 2, we define the model to be used to fit the data. The experimental setup will be detailed in section 3. Section 4 provides a discussion in which alternative models will be compared and also possible uses of this model in the industry. A brief summary and ideas for further research are detailed in section 5.
Numerous formulations of aggregate advertising response models exist in the marketing literature, e.g. little1979aggregate. The model of Nerlove and Arrow nerlove1962optimal extends the Dorfman-Steiner model to cover the situation in which present advertising expenditures affect future demand for products; it is parsimonious and is considered as a standard in the quantitative marketing community. We use it as our starting point.
In this model, advertising expenditures are considered similar in many ways to investments in durable plant and equipment, in the sense that they affect the present and future character of output and, hence, the present and future net revenue of the investing firm. The idea is to define an “advertising stock” called goodwill $A(t)$ which seemingly summarizes the effects of current and past advertising expenditures over demand. Then, the following dynamics is defined for the goodwill
where $u(t)$ is the advertising spending rate (e.g., euros or gross rating points per week), $q$ is a parameter that reflects the advertising quality (an effectiveness coefficient) and $\delta$ is a decay or forgetting rate. The goodwill then increases linearly with the advertisement expenditure but decreases also linearly due to forgetting.
Several extensions and modifications have been proposed to this simple model: it can include a limit for potential costumers vidale1957operations, a non-linear response function to advertise expenditures little1975brandaid, wear-in and wear-out effects of advertising naik1998planning, interactions between different advertising channels bass2007wearout, among other. Still, for most of the tasks, the Nerlove-Arrow model remains as a simple and solid starting point.
Structural time series models or state-space models provide a general formulation that allows a unified treatment of virtually any linear time series model through the general algorithms of the Kalman filter and the associated smoother. Several handbooks durbin2012time, petris2009dynamic, sarkka2013bayesian, west1998bayesian discuss this topic in depth, so we will not develop the corresponding theory. We will however present the most salient features that concern our modelling problem. For further details, the reader may check the aforementioned handbooks.
The state-space formulation of a time series consists of two different equations: the state or evolution equation which determines the dynamics of the state of the system as a first-order Markov process ---\,usually parametrized through state variables\,--- and an observation or measurement equation which links the latent state with the observed state. Both equations are also affected by noise. Denoting by $\mathbf{\theta}_t$ the $m\times1$ state vector describing the inner state of the system, by $G_t$ the $m\times m$ matrix that generates the dynamics, by $H_t$ a $m\times g$ matrix and by $ \epsilon_t $ a $g\times 1$ vector of serially uncorrelated disturbances with mean zero and covariance matrix $W_t$, the system would evolve according to the equation
The states ($\theta_t$) are not generally observable, but are linked to the observation variables $Y_t$ through the observation equation:
where, $Y_t \in \mathbb{R}$ is the observed value at timestep $t$, $F_t$ is the $1\times m$ matrix that links the inner state to the observable, and $V_t \in \mathbb{R}^+$ is the variance of $\epsilon'_t $, the random disturbances of the observations.
The specification of the state-space system is completed by assuming that the initial state vector $\theta_0$ has mean $\mu_0$ and a covariance matrix $\Sigma_0$ and it is uncorrelated with the noise. The problem then consists of estimating the sequence of states $\{\theta_1, \theta_2, \ldots\}$ for a given series of observations $\{y_1, y_2, \ldots\}$ and whichever other structural parameters of the transition and observation matrices. State estimation is readily performed via the Kalman filter; different alternatives however arise when structural parameters are unknown. In the classical setting, these are estimated using maximum likelihood. In the Bayesian approach, the probability distribution about the unknown parameters is updated via Bayes Theorem. If exact computation through conjugate priors is not possible, the probability distributions before each measurement are updated by approximate procedures such as Markov chain Monte Carlo (MCMC) scott2014predicting.
The Bayesian approach offers several advantages compared to classical methods. For instance, it is natural to incorporate external information through the prior distributions. In particular, this will be materialized in Section (ref) where expert information can be incorporated through the spike and slab prior. Another useful advantage is that, due to the Bayesian nature of the model, it is straightforward to obtain predictive intervals through the predictive distribution (see Section (ref)).
The continuous-time Nerlove-Arrow model must be first cast in discrete time so as to formulate our model in state-space. From equation ((ref)), we get:
$$ A_t = (1-\delta) A_{t-1} + q u_{t-1} + \epsilon_t $$ where $A_t$ is the goodwill stock, $u_t$ is the advertising spending rate, $q$ is the effectiveness coefficient and the random disturbance $\epsilon_t$ captures the net effects of the variables that affect the goodwill but cannot be modelled explicitly. This discrete counterpart of Nerlove-Arrow is a distributed-lag structure with geometrically declining weights, i.e., a Koyck model clarke1976econometric, koyck1954distributed. Since in our setting the model includes the effect of $k$ different channels in the goodwill, we shall modify the previous equations according to:
$$ A_t = (1-\delta) A_{t-1} + \sum_{i=1}^k q_i u_{i(t-1)} + \epsilon_t $$
Now, following the notation from ((ref)) and ((ref)), the discrete-time Nerlove-Arrow model in state-space form will read as follows:
The above model is very flexible in the sense that it can be defined modularly, in as much as the different hidden states evolve independently of the others (i.e. the evolution matrix can be cast in block-diagonal form). This greatly simplifies their implementation and allows for simple building-blocks with characteristic behavior. Typical blocks specify trend and seasonal components ---\,which can be helpful to discover additional patterns in the time series--- or explanatory variables that can be added to further reduce the uncertainty in the model and bridge the gap between time series and regression models. Via the superposition principle petris2009dynamic we could include additional blocks in our model:
$$ Y_t = Y_{NA, t} + Y_{R, t} + Y_{T, t} + Y_{S, t} $$ where $Y_{NA, t}$ corresponds to the discretized Nerlove-Arrow equation, defined in ((ref)) and ((ref)); $Y_{R, t}$ is a regression component; containing the effects of external explanatory variables $X_t$; $Y_{T, t}$ is a trend component or a simpler local level component; and $Y_{S, t}$ is a seasonal component.
To take into account the effects of external explanatory variables, a static regression component can be easily incorporated into the model through
$$ Y_{R,t} = X_t \beta + \epsilon_t, $$ where the state $\beta$ is constant over time to favor parsimony.
A spike and slab prior mitchell1988bayesian is used for the static regression component since it can incorporate prior information and also facilitate variable selection. This is specially useful for models with large number of regressors, a typical setting encountered in business scenarios. Let $\gamma$ denote a binary vector that indicates whether the regressors are included in the regression. Specifically, $\gamma_i = 1$ if and only if $\beta_i \neq 0$. The subset of $\beta$ for which $\gamma_i = 1$ will be denoted $\beta_{\gamma}$. Let $\sigma^2_{\epsilon}$ be the residual variance from the regression part. The spike and slab prior george1997approaches can be expressed as $$ p(\beta, \gamma, \sigma^2_{\epsilon}) = p(\beta_{\gamma} | \gamma, \sigma^2_{\epsilon})p(\sigma^2_{\epsilon} | \gamma)p(\gamma). $$ A usual choice for the $\gamma$ prior is a product of Bernoulli distributions:
$$ \gamma \sim \Pi_i \pi_i^{\gamma_i}(1-\pi_i)^{1-\gamma_i}. $$ The manager of the firm may elicit these $\pi_i$ in various ways. A reasonable choice when detailed prior information is unavailable is to set all $\pi_i = \pi$. Then, we may specify an expected number of non-zero coefficients by setting $\pi = k/p$, where $p$ is the total number of regressors. Another possibility is to set $\pi_i = 1$ if the manager believes that the $i-$th regressor is crucial for the model.
Model parameters can be estimated using Markov Chain Monte Carlo simulation, as described in Chapter 4 of petris2009dynamic or scott2014predicting. We follow the same scheme.
Let $\Theta$ be the set of model parameters other than $\beta$ and $\sigma^2_{\epsilon}$. The posterior distribution can be simulated with the following Gibbs sampler:
Repeatedly iterating the above steps gives a sequence of draws $\rho^{(1)}, \rho^{(2)}, \ldots, \rho^{(K)}$ $\sim$ $p(\Theta, \beta, \sigma^2_{\epsilon}, \theta)$. In our experiments, we set $K = 4000$ and discard the first 2000 draws to avoid burn-in issues. \\
In order to sample from the predictive distribution, we follow the usual Bayesian approach summarized by the following predictive equation, in which $y_{1:t}$ denotes the sequence of observed values, and $\bar{y}$ denotes the set of values to the forecast
\[ p(\bar{y} | y_{1:t}) = \int p(\bar{y}| \rho)p(\rho | y_{1:t}) d\rho. \] Thus, it is sufficient to sample from $p(\bar{y}| \rho^{(i)})$, which can be achieved by iterating equations ((ref)) and ((ref)). With these predictive samples $\bar{y}^{(i)}$ we can compute statistics of interest regarding the predictive distribution $p(\bar{y} | y_{1:t})$ such as the mean or variance (MC estimates of $E[ \bar{y} | y_{1:t}]$ and $Var[\bar{y} | y_{1:t}]$, respectively) or quantiles of interest.\\
We can replace the assumption of Gaussian errors with student-$t$ errors in the observation equation, thus leading to the model $$ Y_{t} = F_t \theta_t + \epsilon'_t \qquad \epsilon'_t \sim \mathcal{T}_\nu(0, \tau^2). $$ Typically, in these settings we set $\nu > 1$ to make the variance finite, and this variance parameter can be estimated from data using Empirical Bayes methods, for instance. In this manner, we allow the model to predict occasional larger deviations, which is reasonable in the context of forecasting sales. For instance, a special event not taken into account through the predictor variables may lead to an increase in the sales for that week.
The time series analyzed in this case study contains the total weekly sales of a country-wide franchise of fast food restaurants, Figure (ref), covering the period January 2011 - June 2015, thus comprising 234 observations. The total weekly sales is in fact the aggregated sum from the individual sales of the whole country network of 426 franchises, also allowing a fine-grained study down to the store level, although we shall not carry it out in this paper. Along with the sales figures, the series includes the investment levels $\{u_{it}\}$ in advertising during this period for seven different channels viz. OOH (Out-of-home, i.e. billboards), Radio, TV, Online, \texttt{Search}, \texttt{Press} and \texttt{Cinema}, Figure (ref), $i=1,\ldots,7$.
From such graphs we observe that:
A handful of other predictors $X_t$ which are also known to affect sales will be used in the model, all of them weekly sampled:
Following the notation in Section (ref), we consider three model variants for the particular dataset in increasing order of complexity:
In all cases, only the five principal advertising channels (TV, OOH, ONLINE, SEARCH and RADIO) will be used; the remaining two (CINEMA and \texttt{PRESS}) are sensibly lower both in magnitude and frequency than the others so we can safely disregard them in a first approximation.
As customary in a supervised learning setting with time series data, we perform the following split of our dataset: since it comprises four years of sales, we take the first two years of observations as training set, and the rest as holdout, in which we measure several predictive performance criteria. Before fitting the data, we scale the series to have zero mean and unit variance as this increases MCMC stability. Reported sales forecasts are transformed back to the original scale for easy interpretation. The models were implemented in R using the bsts package scott2016bsts.
It is customary to aim at models achieving good predictive performance. For this reason, we test the predictive performance of our three models using two metrics:
These scores are reported in Table (ref) with sales in million EUR. Note that, unsurprisingly, the models which include external information (RA and RF) achieve better accuracy than the baseline. In addition, we note that predictions are unbiased, since cumulative predictions are extremely close to their observed counterparts. Overall, we found the predictive performance of our models to be successful for a business scenario, as we achieve under 5% relative absolute error using the variants augmented with external information. This is clearly useful for the decision maker who may forecast their weekly sales one week ahead to within a 5% error in the estimation.
Figure (ref) displays the predictive ability of model RF over the hold-out period. Note that the model is sufficiently flexible to adapt to fluctuations such as the peaks at Christmas. Predictive intervals also adjust their width with respect to the time to reflect varying uncertainty, yet in the worst cases they are sufficiently narrow. Further information can be tracked in Figure (ref), where mean standardized residuals are plotted for each model variant. Notice that the residuals for models RA and RF are roughly comparable, being both sensibly smaller than those of the baseline. This means that the simpler Nerlove-Arrow model benefits from the addition of the ambient variables $X_i$, as suggested in our findings from Table (ref).
Having built good predictive models, we inspect them more closely with the aim of performing valuable inferences for our business setting. The average estimated parameter and expected standard deviations of the $q_i$ coefficients for the different advertising channels are displayed in Table (ref). We also show the weights of the ambient variables $X_i$ for the augmented models RA and RF in Table (ref), as well as the probability of a variable being selected in the MCMC simulation for a given model in Figure (ref). Convergence diagnostics of the MCMC scheme are reported in Appendix (ref).
Looking at the ambient variables, the following comments seem in order:
Next, we turn our attention to the investment variables across different advertising channels.
We propose a model which can be used as a decision support system for the manager, helping her in adopting the investment strategy in advertising. The company is interested in maximizing the expected sales for the next period, subject to a budget constraint for the advertising channels and also a risk constraint, i.e., the variance of the predicted sales must be under a certain threshold. This optimization problem based on one-step-ahead forecasts can be formulated as a non linear, but convex, problem that depends on the parameter $\sigma^2$:
where $b_t$ is the total advertising budget for week $t$ and $\sigma$ is a parameter that controls the risk of the sales. We have made explicit the dependence on the regressor variables $x_t$ and advertisement investments $u_t$ in the mean and variance expressions. Solving for different values of $\sigma$, we obtain a continuum of Pareto optimal investment strategies that we can present to the manager, each one representing a different trade-off between risk and expected sales that we can plot in a risk-return diagram. This approach greatly reduces the decision space for the manager.
A possible alternative would be to rewrite the objective function as $$E[ \bar{y}_{t+1} | y_{1:t}, x_{t+1}, u_{t+1}] - \lambda \sqrt{Var[ \bar{y}_{t+1} | y_{1:t}, x_{t+1}, u_{t+1}]}$$ which may be regarded as a lower quantile if the predictive samples $\bar{y}^{(i)}_{t+1}$ are normally distributed. As in the previous approach, different values of $\lambda$ represent different risk-return trade-offs.
If the errors in the observation equation ((ref)) are Normal, the computations for expected predicted sales and the above variance can be done exactly and quickly using conjugacy as in zbMATH06123712. Otherwise, the desired quantities can be computed through Monte Carlo simulation, as in Section (ref).
Note that due to the nature of the state-space model, it is straightforward to extend the previous optimization problem over $k$ timesteps, for $k>1$. The two objectives would be the expected total sales and the total risk over the period $t+1,\dots,t+k$. For normally distributed errors, the predicted sales on the period also follow a Normal distribution, which makes computations specially simple.
The above optimization problem may not need to be solved by exhaustively searching over the space of possible channel investments. In a typical business setting, the manager would consider a discrete set of $S$ investment strategies that are easy to interpret, so she may perform $S$ simulations of the predictive distribution and use the above strategy to discard the Pareto suboptimal strategies.
We have developed a data-driven approach for the management of advertising investments of a firm. First, using the firm's investment levels in advertising, we propose a formulation of the Nerlove-Arrow model via a Bayesian structural time series to predict an economic variable (global sales) which also incorporates information from the external environment (climate, economical situation and special events). The model thus defined offers low predictive errors while maintaining interpretability and can be built in a modular fashion, which offers great flexibility to adapt it to other business scenarios. The model performs variable selection and allows to incorporate prior information via the spike-and-slab prior. It can handle non-gaussian deviations and also provides hints to which of the advertising channels are having positive effects upon sales. The model can be used as a basis for a decision support system by the manager of the firm, helping with the task of allocating ad investments.
Possible extensions of this model include:
The authors acknowledge financial support from the Spanish Ministry of Economy and Competitiveness, through the “Severo Ochoa Programme for Centres of Excellence in R&D” (SEV-2015-0554). V.G. acknowledges support for grant FPU16/05034. The research of D.G-U is supported in part by Spanish MINECO-FEDER Grant MTM2015-65888-C4-3 and MTM2015-72907-EXP. The authors also thank the members of the SPOR-Datalab group at ICMAT for their suggestions and support.