A market model: uncertainty and reachable sets

– Uncertain parameters are always present in models that include human factor. In marketing the uncertain consumer behavior makes it difﬁcult to predict the future events and elaborate good marketing strategies. Sometimes uncertainty is being modeled using stochastic variables. Our approach is quite different. The dynamic market with uncertain parameters is treated using differential inclusions, which permits to determine the corresponding reachable sets. This is not a statistical analysis. We are looking for solutions to the differential inclusions. The purpose of the research is to ﬁnd the way to obtain and visualise the reachable sets, in order to know the limits for the important marketing variables. The modeling method consists in deﬁning the differential inclusion and ﬁnd its solution, using the differential inclusion solver developed by the author. As the result we obtain images of the reachable sets where the main control parameter is the share of investment, being a part of the revenue. As an additional result we also can deﬁne the optimal investment strategy. The conclusion is that the differential inclusion solver can be a useful tool in market model analysis.


Introduction
The market model itself is not the main topic of this article. Though the optimal strategy of the market with respect to the investment is calculated, this is not the main topic, either. The work is focused on the treatment of uncertainty by means of the reachable sets, rather than stochastic modeling. The model we are using (at least the model of the demand) is taken from the literature, with some dynamic properties added.
A lot of demand models that may be used for similar purpose can be found in Lilien and Kotler [1] (linear, non-linear, deterministic, stochastic, static, dynamic models, etc.). What we need is a model that may be adopted into the dynamic market simulation. The demand model we choose is taken from that book. The dynamics is added, to simulate the total accumulated revenue, and the market growth due to the investment. The growth is charged with certain inertia, which implies additional differential equations (see the next section for details).
The search for market models for simulation purpose dates from the early 1960s, and significant publications began to appear in the late 60s. King [2] gives a comparison of iconic, analog and symbolic models. Recall that an iconic model represents reality on a smaller scale, an analog model shows reality by means of maps and diagrams, and a symbolic model uses mathematical expressions to portray reality. Montgomery and Urban [3] discuss the descriptive, predictive and normative models, see Stanovich [4]. By descriptive we mean models that consist largely of diagrams and maps or charts designed to describe a real-world system. Predictive model is used in predictive analytics to create a statistical model of future behavior, and normative model evaluates alternative solutions to answer the question, ''What is going on?'' and suggests what ought to be done or how things should work according to an assumption or standard.
Optimal control in marketing is not a new topic. There are many publications in the field, like, for example Bertsimas and Lo [5]. In that article we can find the application of the Bellman's dynamic programming approach to the problem of the price impact on the dynamic trading on the stock exchange. A tutorial and survey on the relevant technical literature on models of economic growth can be found in Burmeister and Rodney [6]. Feichtinger et al. [7] consider a similar problem for a general market and optimal advertising policy. In that article a detailed formulation of an implementation of the Pontriagin's Maximum Principle is shown. Yuanguo [8] uses the Bellman's Principle of Optimality, and derive the principle of optimality for fuzzy optimal control. This is applied to a portfolio selection model. The applications of differential inclusions (DIs) and reachable set determinations to marketing problems can hardly be found in the known literature. Perhaps the closer topics can be found in my own articles, see Raczynski [17] and Raczynski [9]. However, the problems considered in those articles are quite different from the problem statement of the present research (marketing game with competition and general remarks on DIs in modeling and simulation). Beard [10] presents an application od the DIs to a very specific problem of fertilizers carryover, based on the Liebig's law or the law of the minimum. That law has been developed in agricultural science by Justus von Liebig and states that growth is controlled not by the total amount of resources available, but by the scarcest resource (limiting factor). However, the Beard problem statement is not exactly in the field of marketing and no reachable sets of the model variables are considered.
Read et al. [11] approach the uncertainty in marketing by more heuristic way. Two groups of experts and managers have been asked to think aloud as they make marketing decisions in exactly the same unpredictable situation. The results show significant differences in heuristics used by the two groups. While those without entrepreneurial expertise rely primarily on predictive techniques, expert entrepreneurs tend to invert these. In particular, they use an effectual or non-predictive logic to tackle uncertain market elements and to construct novel markets with committed stakeholders.
Slovik [12] also addresses the problem of market uncertainty. The paper describes the market uncertainty theory comprising the market uncertainty theorem and the notion of heterogeneity of market uncertainty, as well as policy recommendations relevant for rating agencies, financial institutions, and public authorities.
Lien [13] also considers the problem od uncertainty from a cultural and an ethnographic point of view. However, as in the papers of Read and Slovik mentioned above, the problem is stated and described in a textual way only, with no models (logical or mathematical) proposed. This does not mean that their or our approach is better, but means rather that the research of such kinds cannot be compared to each other.
Simos et al. [14] propose a conceptual framework of antecedents and market performance consequences of emergent marketing strategies and test it with a sample of 214 UK enterprises. The results suggest that dimensions of market uncertainty and strategic feedback systems influence the formation of emergent marketing strategy.
Another approach to market uncertainty can be found in Courtney et al. [15]. In that article, four levels of uncertainty are considered: (1) a clear-enough future, (2) alternate futures, (3) a range of futures, and (4) true ambiguity, when multiple dimensions of uncertainty interact and create an environment that is virtually impossible to predict. Compared to a more mathematically rigorous DIs approach, the most similar is the level 3 of uncertainty, where a range of possible future outcomes for marketing data can be known. A heuristic strategies for the different levels of uncertainty are proposed. The general concepts of that work may be useful in the marketing practice, though no formal market model is used.
The article of Brace et al. [16] is an example of one of numerous works on stochastic market models using Heath-Jarrow-Morton (HJM) framework. This, and similar methods are used to analyse forward rates and pricing statistical behavior. It is difficult to compare the DI approach to those methods.
As mentioned also in a comment to the Figure 2 of the present paper, our approach is very different. The main point is that, in marketing and in more general soft system models, not everything can be treated using statistical frameworks.
As stated earlier, the optimization is not our main subject. However, the presented method can also be used to obtain optimal investment strategy. Anyway, to determine reachable sets we use the methods of the optimal control theory.

The model
The demand model we use is a non-linear model of the demand and revenue which parameters are: price, advertising, seasonal index, market overall growth, consumer income, production or acquisition cost and the market elasticities with respect to the price, advertising and consumer income (Lilien and Kotler [1]). The demand model is as follows: where q -demand, q 0 -initial or reference demand, t -the time, s(t) -overall market size, p(t) -price, a(t) -advertising per time unit, y(t) -consumer income, v(t) -seasonal index, e p , e a , e y -market elasticities with respect to price, advertising and consumer income, respectively.
In the model of Lilien and Kotler the market size s(t) is supposed to be equal to (1 + g) t , where g is the market growth rate. We use the variable s instead, to be able to link the market growth to the investment that the company may do in order to expand the market (new installations, infrastructure, etc.).
Normally, the price elasticity is negative (price increase means less demand), and the elasticity for advertising and consumer income are positive. Note that the term v(t)s(t) multiplies the demand, and appears also in the denominator of the advertising impact term. This means that if the market and the seasonal index grow, then we must spend more on advertising to achieve the same effect. We do not consider any storage or warehouse mechanism, so it is supposed that the sales are equal to the demand. The revenue can be calculated as follows: where U is the net revenue, c is the unit cost of the product (production, acquisition cost), a is the advertising, b is the investment factor, p is the price to the public and q is given by the formula (1). The investment b determines what part of the utility e is being invested in the market growth, 0 b 1. The market growth rate is proportional to the product be. It will be one of our control variables. In our equations the market size is relative with initial condition equal to 1. If the company does not invest in the market (b = 0), then the market size has no influence on the model trajectory (remains constant). The final total revenue is the sum (integral) of U over the time interval under consideration. Figure 1 shows the shape of the profit (utility) e as a function of the price and advertising, when the time and other variables are fixed.
It can be seen that the utility function has a maximum that determines the optimal price and advertising at the moment. This is a static case. Now, consider a fixed time interval [0, T ]. To calculate the total revenue we must integrate (2) over the time. The variable b controls the investment, 0 b 1. The market size growth due to the investment. To make the model more realistic, we suppose that the market growth rate follows the investment with certain inertia. The basic equations are as follows: where v ¼ ðp À cÞ qðz; w; y; v; s; tÞ À a; v 0 ¼ initial revenue ds dt ¼ ðx À sÞ=T c market size growth; relative dx dt ¼ g b v=v 0 accumulated investment; relative In the above equations r is the relative accumulated net revenue, s is the relative market size and x is an auxiliary variable. T c is the time constant of the market growth inertia, and g is a constant that defines the impact of the investment in the market expansion. Namely, g tells what is the relative market growth per one currency unit invested. All other variables may change in time. Observe that the additional state variable x has been introduced to manage the inertia.
The final accumulated revenue and the market size strongly depend on the changes of the investment control variable b. So, it is important to be able to see the limits of the revenue and the market size when b changes due to the investment strategy. Note that if b = 1 all the time (maximal investment), then we have no revenue available, and if b = 0 then the market does not grow, which also reduces the available profit.
The aim of the present work is to calculate the reachable sets for all possible investment strategies, as well as to assess the impact of the uncertainty of important market variables.
This should be emphasized that our approach to uncertainty has nothing to do with randomness or stochastic models. If a value of a parameter is uncertain, this means that it can take values from some interval, but this does not mean that it is a random variable. Such a parameter does not have any probability distribution. The changes in its value can be caused by any internal or external agents. For example, in the stock market the information about the actual share price can be obtained by means of observations and predictions, but it also can be a false information introduced intentionally. To obtain the reachable sets, we use the differential inclusion solver, see Raczynski [17,18].

Differential inclusion solver
The DI solver is a software that calculates and displays the solution to a given differential inclusion. The software was developed by S. Raczynski and published in several articles [18]. In Raczynski [9] an application to a marketing differential game is described. However, the topic of that article is different from that presented here. The present research can hardly be compared to other similar published research, because the DI solver had never be used before to marketing models and the shapes of the corresponding reachable sets have not been calculated and presented by anybody (as to the author's knowledge).

Differential inclusions
Let us consider a dynamic model given in the form of an ordinary differential equation (state equation): where x 2 R n is the state vector, t is the time, f is a vectorvalued function and u 2 R m is an external variable, called control in the automatic control theory. Suppose that the value of the control u is restricted so that u(t) 2 C(t), where C(t) & R m is a subset of the R m space. For each fixed x and t, the function f maps the set C (all possible values of u) into a set F & R n . In this way we obtain the following condition.
What we obtained is a differential inclusion (DI). The right-hand side of the equation (4) determines the set F. However, this is merely one of various possible ways to represent the set. In this case it is parameterized by the control u. In general, the DI (5) is not necessarily linked to an ODE model. Let us restrict to a time interval t 2 I = [0, T ]. Any function of time x(t) that obeys the inclusion (5) for almost all t 2 I with a given initial condition, is a trajectory of the differential inclusion. Consider a set Z in I · R n where the graphs of all possible trajectories of the DI that start from a given initial set belong. The set Z is a solution to the DI. The set Z is called reachable set (RS ). Note that a trajectory of a DI is just one possible trajectory, and the solution to a DI with a fixed initial set is given by the reachable set. For more information on differential inclusions consult Aubin and Cellina [19] and Zaremba [20]. In this very natural way we see that the uncertainty in dynamic system modeling leads to differential inclusions as a corresponding mathematical tool. Namely, if any model parameter has uncertain value, it can be treated as a limited control that belongs to a given interval or set. This results in a differential inclusion.
A more extended survey on DI application in modeling and simulation can be found in Raczynski [17]. One of the best publications on DIs in finite-dimensional as well as abstract spaces is the book of Aubin and Cellina [19].

The solver
The DI solver algorithm is not just an extension of any numerical method for the ODEs. We cannot integrate an infinite number of trajectories to obtain the reachable set. One may suppose that RS can be obtained, up to a certain accuracy, integrating a sufficient number of trajectories, each one with different control, generated randomly. Unfortunately, this is not the case and may lead to completely wrong result. The comparison between such random shooting and the DI solver results will be shown in the next section.
One of the properties of the reachable sets is the fact that if a trajectory reaches a point on the boundary of the RS at the final time, then its entire graph must belong to the RS boundary. This fact is well known and used in the optimal control theory (see Pontryagin et al. [21]; Polak [22] and Lee and Markus [23]). Observe that any trajectory that reaches a point on the boundary of the RS is optimal in some sense. Such trajectories can be calculated using several methods, the main one being the Maximum Principle of Pontryagin. If we can calculate a sufficient number of trajectories that scan the RS boundary, then we can see its shape. This can be done by some kind of random shooting over the RS boundary. Such shooting has nothing to do with a simple or primitive random shooting mentioned earlier, when the trajectories are generated randomly inside the RS.
In few words, the DI solver works as follows. The user provides the DI in the form of an equivalent control system. To do it he/she must parameterize the right-hand size (the set F) using a m-dimensional auxiliary control variable u. The DI solver automatically generates the equations of the so-called conjugate vector p and integrates a set of trajectories, each of them belonging to the boundary of the RS. To achieve this, over each trajectory the Hamiltonian H(x, p, u, t) is maximized. This procedure is similar to that used in dynamic optimization. In the optimal control problem the main difficulty consists in the boundary conditions for the state and conjugate vectors. For the state vector we have the initial condition given, and for the conjugate vector only the final conditions (at the end of the trajectory) are known, given by the transversality conditions. This means that the optimal control algorithm must resolve the corresponding two-point-boundary value problem. In the case of a DI we are in a better situation. There are no objective function and no transversality conditions. As the consequence, for the vector p we can define the final as well as the initial conditions. Anyway we obtain a trajectory which graph belongs to the RS boundary. With the initial conditions for p we integrate the trajectory only once, going forward. The only problem is how to generate these initial conditions in order to scan the RS boundary with nearly uniform density. The algorithm is quite simple: the initial conditions for p are generated randomly, due to a density function that is being automatically modified, covering with more density points that correspond to trajectories that fall into low density spots at the RS boundary. Trajectories that are very close to each other are not stored (storing only one from each eventual cluster). As the result we obtain a set of trajectories covering the RS boundary that can be observed in graphical form and processed.
Application to the market model Note that our model (3) has the form of a control system, and defines the corresponding DI. In the investment strategy problem, the control variable is b(t), the share of investment we apply out of the total revenue (Experiment 1, below). The uncertainty of the other model variables can be treated in a similar way, letting the variable change within a given interval. Below are presented the results of such analysis for the uncertainty in the values of e p , -the market elasticity with respect to the price (Experiment 2).

Experiment 1
In this case our aim is to obtain the reachable set in the revenue/market size plane for all possible strategies of investment, with final time fixed. As an additional result, the DI solver provides the control function (investment as a function of time) for any point on the boundary of the reachable set. So, we can obtain the optimal strategy that maximizes the total revenue or the final market size.
The model is given by the equations (3), where b is the control variable, that can change within 0 and 1. The model parameters are as follows: Initial sales = 70 000 Initial (reference) advertisement = 10 000 per time unit (a day) Product unitary cost = 0.8 Initial relative market size = 1 Unit price = 1 Unit initial relative price = 1.2 Time constant for investment inertia T c = 30 days g = 0.006 Market elasticity for the product price = À2 Market elasticity for advertisement impact = 0.6 Final simulation time = 365 days Note that the initial relative price is set equal to 1.2 and not 1. This is to make the demand change when the price elasticity changes. In this experiment both seasonal index and the consumer income are functions of time. The relative consumer income is equal to 1 except the interval (200, 214) where it is equal to 2.5. The seasonal index is equal to 3 within the interval (100, 128), equal to 0.5 in (160, 167) and equal to 1 elsewhere. Figure 2 shows the image of the reachable set for the revenue and market size at the simulation final time, with the investment control b included in [0, 1]. It should be noted that in our case the dimensionality of the state space is equal to 3. So, to show the complete reachable set we should display it in the 4-dimensional time-state space, which is quite difficult. If we limit the displayed state space components to market size-revenue, then the reachable set boundary is not always seen just one contour, because for each time instant the 2-dimensional image is a projection of a 3-dimensional reachable region. Fortunately, in this case the boundary is clearly seen, perhaps because the control vector dimensionality is equal to 1. On the same figure we also can see a set of reachable points obtained by simple random shooting (10 000 trajectories integrated). It is clear that the simple shooting provides a wrong assessment of the reachable set.
The result of Figure 2 shows the advantage of the DI approach compared to any other method which uses models with random variables to treat uncertainty (stochastic approach). Whatever the distribution of the uncertain parameter would be, the estimate of the reachable set will result to be wrong. Moreover, recall that uncertain parameters are not necessarily random ones. An uncertain parameter may have no probability distribution at all, and may have quite deterministic value, like, for example a false information generated intentionally (this occurs frequently in the stock market). Figure 3 depicts the 3D image of the reachable set in the time-state space.
After the set has been displayed, the user can select any point of its boundary and see the corresponding control. Selecting the point with maximal revenue we obtain a control of ''bang-bang'' type, equal to 1 for time between 0 and 97 and equal to zero between 97 and 365 days. The changes in the corresponding accumulated net revenue are shown on Figure 4.

Experiment 2
Now, suppose that the value of the market elasticity with respect to the price e p is uncertain and may change up to plus minus 20% of its original value. Using model parameters as in Experiment 1, with variable investment control in [0, 1] and uncertain (variable) e p , we obtain a larger reachable set. In this case the set is more complicated (2D projection of a 3D region).
This results in a cloud of reachable points and not in a simple contour. A post-processing of this cloud provides an assessment of the reachable set, shown as the gray region on Figure 5. To compare with Experiment 1, the boundary points of the previous reachable set are also shown in the same figure.
Observe both the net revenue and the market size may reach greater values than that obtained with fixed elasticity e p .
In general, it can be seen that our market needs the optimized investment strategy to reach maximal revenue. Note that we have fixed final simulation time, as if the company sells   during the given time interval and then disappears. Experiments with moving time horizon can be done, but may be more complicated and time-consuming.

Experiment 3: Uncertainty and reachable set vs. sensitivity
Our approach to uncertainty is based on model reachable sets. As stated before, this is a deterministic approach. We do not deal with any probability distribution for the uncertain parameters, using rather the intervals where the parameters may belong. The influence of the parameters is being treated dynamically, i.e. we allow that the parameters are functions of time, changing over each model trajectory. The reachable sets obtained in this way show the influence of the uncertain parameters, somewhat similar to the sensitivity analysis. Let us recall that the basic sensitivity analysis shows the possible changes in the model state due to the changes of the initial conditions and the (constant) model parameters. Our problem statement and that of sensitivity methods are quite different to each other, so it is difficult to compare these methods. We do not pretend to prove if the DI applications are better or worse from the results of sensitivity methods. It is just interesting to compare them.
The sensitivity analysis is not the topic of this paper, so we will not provide any review of the (huge) literature in the field. For a comprehensive review see Hamby [24]. The common way to conduct sensitivity analysis is to assign probability density functions to each model parameter and assess the resulting outcome variance with respect to the corresponding parameters. Though the advanced methods is sensitivity research are quite sophisticated and perhaps more accurate, we will use in this comparison a basic sensitivity concepts.
One of the commonly known methods is the general error propagation formula, i.e.
where Y is the model output (the outcome being analyzed), V(Y) is the variance of Y, the model parameters are X 1 ,. . ., X n , and V(X i ) is the variance of the parameter X i . In our case Y may be the value of any of the model state variables at the end point of the simulated trajectory. We do not have any algebraic expression for Y, which value is the result of model simulation. However, as the calculation of a single trajectory take only a fraction of a second, the partial derivatives can be easily calculated by numerical differentiation.  Figure 6. Bold rectangle shows the region of plus minus standard deviation for the outcome. It can be seen that this region is little informative, and obviously very different from the reachable set (gray area). Formula (6) has been obtained using model linearization. As our model is non-linear, a better, direct way to see the influence of the parameters is simply integrating model equations, changing the parameters, each of them being constant over the trajectory. The result is also shown on Figure 6 as a dotted area.   Observe that the true reachable set is greater than the region obtained this way. For more complex models, in particular for oscillating systems, this difference is greater.

Conclusions
Calculating reachable sets is something more than obtaining the optimal control. Though we can get the optimal control strategy as an additional result, the possibility of viewing the shape of reachable sets provides much more information about the model behavior. The same simulation program can be used for many more experiments, looking for reachable sets with variable advertising and the uncertainty of any one of the model parameters. Adding more inertia to the model will result with higher state vector dimensionality and more complicated images of the reachable sets. For such sets the problem is not only the set determination, but also the way to display it. As the computer screen is (still) 2-dimensional, the good effects can be obtained while displaying the cloud of points in movement, rotating around one of the state space axes. Unfortunately, such moving displays can hardly be shown as figures of a printed article. The comparison of the DI approach with sensitivity analysis is interesting, but it should be remembered that the problem statement in both approaches is quite different. Most of the sensitivity methods are designed to find the variance propagation and statistical properties of the model output, while the DIs and the reachable sets are deterministic. The selection of an appropriate methodology always depend on which is the purpose of the research.