The Markets of Green Cars of Three Countries: Analysis Using Lotka–Volterra and Bertalan ﬀ y–Pütter Models

: Did the diesel scandal of 2015 a ﬀ ect the market for cars? We consider this question in relation to Germany, Austria, and Switzerland. Starting with historical registration data of cars with di ﬀ erent drivetrain technologies, we considered each technology in isolation and ﬁtted a ﬁve-parameter Bertalan ﬀ y–Pütter (BP) growth model to the stocks of cars. We used this model as it generalizes several well-known three-parameter models, which are distinguished by their exponent pair, e.g., Brody model BP (0, 1), West model BP (0.75, 1), and logistic growth BP (1, 2). We then used these models to derive a Lotka–Volterra (LV) model for the co-evolution of the (annual) market shares of the di ﬀ erent drivetrain technologies. We augmented this model by a consideration of model uncertainty and found that initially all technologies were in a state of competition, except for Austria, which changed in 2015 to a predator–prey situation with diesel as the sole prey. This analysis of model uncertainty compared the best-ﬁtting growth curve with the growth trajectories of other likely (Akaike weight 5% or higher) models of BP type. We conclude with remarks about open innovation.


Background
The use of conventional cars in cities has been criticized for hazardous emissions of nitrogen oxide (NOx) and particulate matter [1]. These emissions are known to trigger asthma, bronchitis, and other respiratory and cardiovascular diseases [2,3]. In response, authorities worldwide have made emission standards more stringent in order to pressure the automotive industry towards developing more environmentally friendly cars.
This legislation raised interest in environmentally friendly alternatives, such as hybrid cars, fuel-cell cars, and battery-driven electric cars [4]. For instance, China has formulated a strategic plan to promote market penetration of battery-powered electric vehicles, aiming at resolving the massive pollution (haze) of its mega-cities as well as reducing dependence on imported oil [5]. In Europe, the green cars initiative developed a roadmap toward electrification, pointing out that the major differences between conventional and green technologies would be the aspects of energy and resource security, climate change, public health, freedom of mobility, and economic growth [6]. The European Commission, too, noted that transport sector electrification would be essential "for meeting the European Union goals of decarbonization and energy security, as it accounts for 25% of all CO 2 emissions in Europe" [7].

Outline of the Approach
The progress in the analytic methods for market diffusion is illustrated by the proliferation of models, but there is no consensus on the best of these approaches [22]. In this paper, we aim at modeling the growth of the stock of cars with different drivetrain technologies (Table 1). For this type of problem, trend models, such as logistic growth, are common [5,23]. We improve this by using the Bertalanffy-Pütter (BP) growth model, which generalizes the common trend models and therefore achieves a better fit. In the initial step, we model the growth curves of each drivetrain technology independently of the other technologies by using the BP model to describe its stock (y) of cars over time (t). The growth functions y(t) of the BP model are solutions of the differential Equation (1) of Pütter [24]. The parameters of the model are determined from fitting the growth curve to data (details below and in Table 2): We then ask if perhaps a simpler model might suffice, using three instead of five parameters; examples are the logistic and the Gompertz growth models.
Trend models have the disadvantage that they investigate single products in isolation. We therefore continue with a system dynamic approach [28] that uses a generalized Lotka-Volterra (LV) system of differential Equation (2) to model the dynamics of the market shares of the competing drivetrain technologies. We use this model to unveil otherwise hidden changes in the characteristics of the market dynamics of the three countries. Thereby, one of the technologies (index i = 0) is the outside, which supposedly does not affect the market dynamics. For the other technologies, the market shares z i (t) interact as in Equation (2). We will use the BP model to find suitable coefficients g i (t) for the LV equations (details below).

Data
For Table 1, we retrieved the annual registrations (sales figures of new cars) from publicly available online resources. We distinguished gasoline-powered cars, diesel cars, electric cars, hybrid cars, and cars powered by natural gas or liquid gas (liquified petroleum gas). To make the data comparable, we combined certain data (e.g., hybrid cars using gasoline or using diesel) and we ignored rare technologies (e.g., Switzerland counted vehicles without a motor). Note that, for Switzerland, the data for hybrid cars were not available for the first year.
As growth models in general are not suitable in situations where growth and decay both occur, we modeled the stock of cars. We therefore used Table 1 to compute the accumulated registration data for the mentioned technologies. If x 1 , x 2 , x 3 , . . . are the registrations from successive years t 1 , t 2 , t 3 , . . . (Table 1), then the stock is y 1 = x 1 , y 2 = y 1 + x 2 , . . . , y k = y k−1 + x k (sum of new sales from the first year until present). We thereby also achieve a smoothening of the data. For the modeling, we represent the years by t 1 = 0, t 2 = 1, t 3 = 2, . . . (except for hybrid CH: t 1 = 1, as the other CH time series started one year earlier).

BP Model
We use the BP growth model, differential Equation (1), to describe the stock y(t) of cars over time. Equation (1) can be solved analytically, though, in general, not by means of elementary functions [30]. The parameters of the model (Table 1) are the exponent pair a < b, non-negative constants p and q, and the initial value y(0) = c > 0, where, for each country, t = 0 is the first year of its data (t = 1 for Swiss hybrid cars). They are determined from fitting the growth curve to data; see Table 2. As we optimize also for the initial values c, the best-fit values for c in Table 2 will slightly differ from the observed initial registration data in Table 1.

Other Growth Models
In the literature, there are several other five-parameter growth models, such as the model of Bass [31] for market diffusion, which has also been used to model the diffusion of electric cars [32]. We decided to use the BP model as it is known to be flexible [33] and as it includes as special cases several simple three-parameter models which have been used previously in the analysis and forecasting of technology diffusion and business trends [34,35]. Specifically, the modeling of the automotive market has often used logistic growth [5,23] or the Gompertz model [36][37][38].
The Gompertz [45] model is the limit case BP(1, 1) with a different differential equation, where b converges to a = 1 from above [46]. When we compare the Gompertz model with other BP models, we use BP(1, 1.01) as a proxy for the Gompertz model. Figure 1. Exponent pairs as proxy for models. Note: search grid (yellow), exponent pairs with Akaike weight 5% or higher (red), best-fit exponent pair (black), and exponent pairs of selected named models (blue). Plot using Mathematica 12.1 based on the stock of electric cars in Austria (accumulated registration data computed from Table 1).

Data
For Table 1, we retrieved the annual registrations (sales figures of new cars) from publicly available online resources. We distinguished gasoline-powered cars, diesel cars, electric cars, hybrid cars, and cars powered by natural gas or liquid gas (liquified petroleum gas). To make the data comparable, we combined certain data (e.g., hybrid cars using gasoline or using diesel) and we ignored rare technologies (e.g., Switzerland counted vehicles without a motor). Note that, for Switzerland, the data for hybrid cars were not available for the first year.
As growth models in general are not suitable in situations where growth and decay both occur, we modeled the stock of cars. We therefore used Table 1 to compute the accumulated registration data for the mentioned technologies. If x1, x2, x3, …are the registrations from successive years t1, t2, t3, ... (Table 1), then the stock is y1 = x1, y2 = y1 + x2, …, yk = yk−1 + xk (sum of new sales from the first year until present). We thereby also achieve a smoothening of the data. For the modeling, we represent the years by t1 = 0, t2 = 1, t3 = 2, ... (except for hybrid CH: t1 = 1, as the other CH time series started one year earlier).

BP Model
We use the BP growth model, differential Equation (1), to describe the stock y(t) of cars over time. Equation (1) can be solved analytically, though, in general, not by means of elementary functions [30]. The parameters of the model (Table 1) are the exponent pair a < b, non-negative constants p and q, and the initial value y(0) = c > 0, where, for each country, t = 0 is the first year of its data (t = 1 for Swiss hybrid cars). They are determined from fitting the growth curve to data; see Table  2. As we optimize also for the initial values c, the best-fit values for c in Table 2 will slightly differ from the observed initial registration data in Table 1.

Other Growth Models
In the literature, there are several other five-parameter growth models, such as the model of Bass [31] for market diffusion, which has also been used to model the diffusion of electric cars [32]. We decided to use the BP model as it is known to be flexible [33] and as it includes as special cases several simple three-parameter models which have been used previously in the analysis and forecasting of Figure 1. Exponent pairs as proxy for models. Note: search grid (yellow), exponent pairs with Akaike weight 5% or higher (red), best-fit exponent pair (black), and exponent pairs of selected named models (blue). Plot using Mathematica 12.1 based on the stock of electric cars in Austria (accumulated registration data computed from Table 1).

Sales and Market Shares
Starting with the BP model curve y i (t) fitted to the stock of cars using a certain drivetrain technology (i), we use the derivative y i as the model curve for the current annual sales . For simplicity, we replace ζ by a fixed time shift 0 ≤ ts i ≤ 1 for all k ≥ 1; i.e., we model sales x k+1 by the function y i (k − ts i ).
To obtain the time shift ts i , we use the method of least squares, minimizing the sum of squared errors between x k+1 and y i (k − ts i ), where y i is the best-fit BP growth function; see Table 4. The market share s i (t) at time t is then given by Equation (3), where the total current sales are the market volume m(t): and The computations related to the market shares for Germany use the combined data for natural and liquid gas cars. This ensures that the market dynamics of different countries are compared using the same number of comparable technologies

LV Model
Lotka [47] and Volterra [48] introduced a system of differential equations to model the dynamics of competing species. Despite its simplicity, this model displays a multitude of possible patterns of co-existence [49] that made it an attractive model of economic competition [50,51]. Here, we use the Lotka-Volterra (LV) system (2) with variable coefficients, as, for this system, there is an explicit solution [52,53]. Variable coefficients are used as they represent situations where companies may actively influence the competitive roles of their products (e.g., PR and marketing, technological advances).
One of the technologies (index i = 0) is selected as the outside. Generally, this is the technology with the lowest market share. For the other technologies, we use the market shares s i (t) from above to define certain functions f i and g i ; see Equation (4). For i ≥ 1, the market shares are then solutions s i (t) = z i (t) of the LV system (2) with these coefficients g i . The market share z 0 (t) is the difference to 100% of the sum of the z i (t), i ≥ 1 (remainder of the market); z 0 = s 0 .
This model has an economic meaning [54]: f i are the (assumed) utilities of consumers for different drive technologies (the utility of the outside good i = 0 is 0) and, if consumers maximize their utility, then Equation (2) describes the market shares. Amongst economic applications are, e.g., studies of the competition between tourism destinations [55] or between ports [56,57]. We use this model to identify the competitive market characteristics [52]. If the two coefficients g i (t) and g j (t) are positive, then, at time t, the two technologies (i, j) are in a state of competition. If all coefficients g i (t) of Equation (2) are positive, this is a state of pairwise competition. This means that (at time t) more sales of cars using any other technology reduce the growth of the own market share. In the absence of other technologies, the own market share would grow. If one of the two coefficients g i (t) or g j (t) is positive and the other is negative, then the two technologies are in a predator-prey situation. The technology with the negative coefficient is the prey, the other one the predator, whereby more prey enhances the growth of the predator and more predators inhibit the growth of prey. (Other than in the classical ecological predator-prey model, in this model, the market shares of predators/prey would grow/decrease also in the absence of prey/predators, whence the process is not circular.) If both coefficients g i (t) and g j (t) are negative, then the two technologies are in a situation of mutualism (symbiosis), where the presence of the other technology enhances the own growth. The exceptional cases, where one or both coefficients vanish, are characterized as commensalism, amensalism, and neutralism. These cases were observed for LV models of innovation diffusion [58].

Calibration
We use various measures for the goodness of fit. They are all related to the sum SSE of squared errors: when fitting the BP growth function to the stock of cars, we seek parameters that minimize SSE. Thereby, if y(t) is a solution of Equation (1), using certain exponents a < b and parameters p, q, c, if y k are the stock data at time t k , and if n is the number of data, then SSE is defined by Equation (5).
In Equation (6), RMSE is the root-mean-squared-error and NRMSE is the normalized RMSE (as percent of the mean). R-squared (R 2 ) of Equation (7) is the coefficient of determination that compares the model fit with the fit by the trivial constant model.

Data Fitting
Standard optimization tools (e.g., NonlinearModelFit of Mathematica) may not always find best-fit parameters for model (1). The literature has reported problems already for the simpler Richards model [59]. This paper uses a grid-search strategy for data fitting. It identifies the best-fit parameters for 4 × 10 4 three-parameter models BP (a, b) using exponent pairs (a, b) on a certain grid, the yellow area in Figure 1. Figure 2 plots these growth curves (they cover the yellow area) together with the best-fit curve corresponding to the best-fit exponent pair (black).  Table 2, band of growth curves with an Akaike weight of 5% or higher (red), and band of the best-fit growth curves using any exponent pair from the search grid (yellow); plot using Mathematica 12.1.
The search starts with a grid (step size 0.01 in both directions) of exponent pairs ( Figure 1). For each grid point (a, b), we seek the best-fit parameters (p, q, c) for the three-parameter model BP (a, b); the so minimized sum of squared errors is SSEopt (a, b). This optimization uses a custom-made variant of the method of simulated annealing [60,61]. The outcome of optimization is exported into a spreadsheet (its columns list a, b, c, p, q, and SSE). The best-fit parameters (amin, bmin, cmin, pmin, qmin) of the grid are identified from the row with the lowest value of SSE. This defines the overall lowest SSE of (1) as SSEmin = SSEopt (amin, bmin). If the best-fit exponent pair is on the upper or left boundary of the grid, then we enlarge the grid near this boundary and repeat the search, until an optimal inner point is found. This optimizes the exponents amin, bmin with an accuracy value of 0.01 (for the other parameters, we use a higher accuracy value.)

Model Uncertainty
For each type of technology and each country, there is a different best-fit model for the growth of the stock of cars. We now ask if a simpler modeling approach is feasible: is there one exponent pair (a, b), so that BP (a, b) has a good fit to all growth data? For instance, can all growth curves be modeled by the often-used logistic model? We use two approaches to address this question. One measures goodness of fit in terms of R 2 , the coefficient of determination. It asks which R 2 values were obtained for the data using a specific model.
The second approach focuses on the models that are likely to be true. In Equation (8), AIC is the Akaike [29] information criterion; it penalizes models that use more parameters (overfitting), whereby the number p of parameters also counts SSE (e.g., p = 4 for the logistic model). Given a finite set of models, the Akaike weight prob compares a given model with the overall best-fit model (it has the lowest AIC: AICmin). Assuming that either the best-fit model or the given model is true, then prob is the probability (maximum: 50%) that the given model is true despite its poorer fit [62]. As in Figure  1 (red area), we then identify all models with an Akaike weight of 5% or higher.  Table 2, band of growth curves with an Akaike weight of 5% or higher (red), and band of the best-fit growth curves using any exponent pair from the search grid (yellow); plot using Mathematica 12.1.
The search starts with a grid (step size 0.01 in both directions) of exponent pairs (Figure 1). For each grid point (a, b), we seek the best-fit parameters (p, q, c) for the three-parameter model BP (a, b); the so minimized sum of squared errors is SSE opt (a, b). This optimization uses a custom-made variant of the method of simulated annealing [60,61]. The outcome of optimization is exported into a spreadsheet (its columns list a, b, c, p, q, and SSE). The best-fit parameters (a min , b min , c min , p min , q min ) of the grid are identified from the row with the lowest value of SSE. This defines the overall lowest SSE of (1) as SSE min = SSE opt (a min , b min ). If the best-fit exponent pair is on the upper or left boundary of the grid, then we enlarge the grid near this boundary and repeat the search, until an optimal inner point is found. This optimizes the exponents a min , b min with an accuracy value of 0.01 (for the other parameters, we use a higher accuracy value.)

Model Uncertainty
For each type of technology and each country, there is a different best-fit model for the growth of the stock of cars. We now ask if a simpler modeling approach is feasible: is there one exponent pair (a, b), so that BP (a, b) has a good fit to all growth data? For instance, can all growth curves be modeled by the often-used logistic model? We use two approaches to address this question. One measures goodness of fit in terms of R 2 , the coefficient of determination. It asks which R 2 values were obtained for the data using a specific model.
The second approach focuses on the models that are likely to be true. In Equation (8), AIC is the Akaike [29] information criterion; it penalizes models that use more parameters (overfitting), whereby the number p of parameters also counts SSE (e.g., p = 4 for the logistic model). Given a finite set of models, the Akaike weight prob compares a given model with the overall best-fit model (it has the lowest AIC: AIC min ). Assuming that either the best-fit model or the given model is true, then prob is the probability (maximum: 50%) that the given model is true despite its poorer fit [62]. As in Figure 1 (red area), we then identify all models with an Akaike weight of 5% or higher.

Simulation
Finally, we explore the extent to which the growth trajectories depend on the choice of growth model (model uncertainty). To this end, we use a Monte Carlo simulation to identify the growth patterns of the likely growth curves (red area in Figure 2). To avoid a "junk in-junk out" situation from a simulation with generally false models, we exclude models from the simulation that are unlikely to be true. Thus, we select at random 5000 exponent pairs (a, b) of the grid (red area in Figure 1) so that the model BP(a, b) has an Akaike weight (probability to be true) of prob > 5% when compared to the best-fit model BP(a min , b min ). We therefore do not penalize the best-fit model for the optimization of the two exponents, as the optimization was merely a comparison of a given finite set of three-parameter models defined by a search grid. (Thus, in formula (8), we use p = 4 to count the parameters c, p, q, and SSE; n is the number of time-points.) The upper and the lower function values of the 5000 best-fit growth curves then define the band of the growth scenarios that can be obtained from the likely models (Figures 2 and 3

Simulation
Finally, we explore the extent to which the growth trajectories depend on the choice of growth model (model uncertainty). To this end, we use a Monte Carlo simulation to identify the growth patterns of the likely growth curves (red area in Figure 2). To avoid a "junk in-junk out" situation from a simulation with generally false models, we exclude models from the simulation that are unlikely to be true. Thus, we select at random 5000 exponent pairs (a, b) of the grid (red area in Figure  1) so that the model BP(a, b) has an Akaike weight (probability to be true) of prob > 5% when compared to the best-fit model BP(amin, bmin). We therefore do not penalize the best-fit model for the optimization of the two exponents, as the optimization was merely a comparison of a given finite set of threeparameter models defined by a search grid. (Thus, in formula (8), we use p = 4 to count the parameters c, p, q, and SSE; n is the number of time-points.) The upper and the lower function values of the 5000 best-fit growth curves then define the band of the growth scenarios that can be obtained from the likely models (Figures 2 and 3, red area). We apply the simulation approach to identify the most likely scenarios for the competition between the different technologies. To this end, for each technology, we select at random a growth curve with Akaike weight of at least 5% and we combine these growth curves to define the market shares and the corresponding LV equations. We then classify for each simulated dynamic model its pattern of competition (scenario) at a given moment in time. The outcome indicates the percentage of the 5000 simulations for which a certain scenario is observed. We are particularly interested in three scenarios: pair-wise competition of all technologies (meaning business as usual), diesel as the sole prey (perhaps a consequence of the diesel scandal), and diesel together with petrol as the sole prey (a preference of consumers for green technologies). In addition, we inform about the relative weights of the scenarios. Therefore, the weight of a simulation is the product of the Akaike weights of the selected growth functions for the different technologies. The total weight is the sum of the weights of all 5000 simulations. We apply the simulation approach to identify the most likely scenarios for the competition between the different technologies. To this end, for each technology, we select at random a growth curve with Akaike weight of at least 5% and we combine these growth curves to define the market shares and the corresponding LV equations. We then classify for each simulated dynamic model its pattern of competition (scenario) at a given moment in time. The outcome indicates the percentage of the 5000 simulations for which a certain scenario is observed. We are particularly interested in three scenarios: pair-wise competition of all technologies (meaning business as usual), diesel as the sole prey (perhaps a consequence of the diesel scandal), and diesel together with petrol as the sole prey (a preference of consumers for green technologies). In addition, we inform about the relative weights of the scenarios. Therefore, the weight of a simulation is the product of the Akaike weights of the selected growth functions for the different technologies. The total weight is the sum of the weights of all 5000 simulations.

Graphical Review of Methodology
Figures 1-5 explain our methodological approach for a concrete example, the electric car data from Austria ( Table 1). The grid search has reduced optimization to the search for the best-fitting model amongst those with exponent pairs in the search grid (yellow area in Figure 1). The blue dots in Figure 1 are the exponent pairs of certain named three-parameter growth models. The black dot is the exponent pair (a min , b min ) of the best-fit model. None of the exponent pairs of the named models have been close to this optimal exponent pair. The red area of Figure 1 collects 22,442 exponent pairs of BP models with an Akaike weight of 5% or higher when compared to the best-fit model. The frayed appearance of the red area reflects the random character of simulated annealing; it leads to slight fluctuations of the accuracy for different grid points.

Graphical Review of Methodology
Figures 1-5 explain our methodological approach for a concrete example, the electric car data from Austria ( Table 1). The grid search has reduced optimization to the search for the best-fitting model amongst those with exponent pairs in the search grid (yellow area in Figure 1). The blue dots in Figure 1 are the exponent pairs of certain named three-parameter growth models. The black dot is the exponent pair (amin, bmin) of the best-fit model. None of the exponent pairs of the named models have been close to this optimal exponent pair. The red area of Figure 1 collects 22,442 exponent pairs of BP models with an Akaike weight of 5% or higher when compared to the best-fit model. The frayed appearance of the red area reflects the random character of simulated annealing; it leads to slight fluctuations of the accuracy for different grid points.    Figure 2 plots the stock (accumulated registration data) of electric cars (AT) during the years 2011-2019 (blue) and the best-fit BP growth curve to these data (black). It corresponds to the black exponent pair and uses the best-fit exponents from Table 2. The yellow and red bands in Figure 2 summarize the outcomes of Monte Carlo simulations: it displays, for discrete points in time (t = 0, 1, . . . , 10), the minimum and maximum of the function values of the 5000 randomly selected growth functions whose exponent pairs come from the yellow and red areas of Figure 1, respectively. The figure connects these points by lines and the intermediate area is shaded accordingly. The yellow area is relatively broad, as it includes also poorly fitting models. However, as each of the randomly selected model curves is a best-fit curve subject to the choice of the exponent pair, the yellow area is not too remote from the data-points. The red band (growth curves of other likely models) is close to the data.     Table 1 for electric cars (AT) with the derivatives of the best-fit curve and the curves from the random sample with red exponent points (Figure 1). Table 4 lists the used time shifts. For the simulations, we used the same time shift ts as for the best-fit model. For the annual sales data (x 2 , x 3 , . . . , x 10 of Table 1), the fluctuations around the model curve y (t-ts) are much higher than the fluctuations of the accumulated data y n around the BP growth curve y(t). However, as the red band (sales data from other likely models) is wide, high fluctuations are to be expected. Figure 4 compares the market shares for electric cars (AT) with Formula (2), the estimate for the market shares using the best-fit model. As the market share cannot be computed in isolation, this computation uses the best-fitting growth curves for all technologies. The computation of the red band of the most likely market shares uses 5000 simulations, where, for each technology, one best-fitting BP-growth curve is selected at random. More specifically, for each technology, one exponent pair (a, b) in the red region (of the analogue to Figure 1) is selected at random and the growth curve from the model BP(a, b) with the best fit to the data of the technology is used. Figure 5 plots the consumers' utility for electric cars (AT) based on the best-fit models for the market shares with gas cars as the outside goods. It compares this curve to the range of the other likely utilities. A common observation for all best-fit utility functions for all countries was their approximately parabolic shape.

Best-Fit Exponent Pairs
Several authors have modeled the growth of the car market by the logistic model [5,23]. This model achieved an acceptable fit for all 16 data (R 2 > 95.21%). The BP models of Table 2 improved the fit significantly. None of the exponent pairs of a named model were close to a best-fit exponent of Table 2; an exception was the Gompertz model (close to electric car GE).
From Table 2, the following pattern emerged that seems to relate the technology of a car to the location of the exponent pair of the best-fitting BP model. The first exponent a was small for gas cars and the difference b-a was low for electric and hybrid cars. Indeed, for three gas car data (natural gas in AT and CH, liquid gas in GE), the first exponent was a = 0 (the growth curve was not sigmoidal). For the fourth gas car data (natural gas in GE), the first exponent was a = 0.01. For the three hybrid car data and for two electric car data, the best-fit exponents were close to the diagonal (b-a = 0.01). For the third electric car data (AT), the difference between the exponents was 0.12. Furthermore, b-a was high for diesel cars; the maximal differences of the exponents (b-a > 2.9) were realized for two diesel data (AT, CH).
The original motivation for the BP model was biological [63]: the exponent pairs characterize certain biophysical situations about the growth of animals [42,43]. While this suggests that the hypothesis that animals of the same species have similar best-fit exponent pairs to describe their growth, in our previous work, we did not find evidence to support this hypothesis. We therefore were surprised by the result of this paper that the best-fit exponents for modeling the growth of the stock of cars from different countries allowed us to distinguish certain drivetrain technologies: gas cars from all three countries were characterized by a small first exponent. Green cars (electric or hybrid) were characterized by a small difference between the exponents.

Exponent Pairs Leading to Reasonable Fits
Seeking one BP model with a reasonable fit to all data in terms of R-squared, the best choice was BP(0.7, 0.88), with R 2 > 96.93% for all 16 data; this was a slight improvement over the logistic model. Much better fits were achieved (Table 5) if, for each drivetrain technology, a common model was sought. The reason was antagonism between gas-powered cars and electric or hybrid cars. (To achieve a better fit, for the latter, exponent pairs close to the line a = 0 were needed, while, for the former, the exponent pairs with a better fit were more remote from this line. We observed this antagonism for each country.) The exponent pairs of these models display a similar pattern with respect to the drivetrain technology, as we observed previously for the best-fit exponent pairs. Furthermore, while for nine of the 16 data, an excellent fit of R 2 > 99.8% could be achieved (Table 2), no single three-parameter BP model achieved this good fit for all nine data. The best outcome was achieved by BP(1.19, 1.2), with R 2 > 99.8% for seven data. Furthermore, for different technologies, certain named trend models provided fits that were close to the best achievable fits (Table 5): the Brody model was suitable for gas cars, the Gompertz model (logistic growth was suitable, too) for electric and hybrid cars (this confirms the finding of [36]), and the West model (the Bertalanffy model and logistic growth were suitable, too) for conventional cars (diesel, gasoline). For these models, R-squared was close to the largest R-squared that could be achieved by a single model for all datasets of these classes. This observation was corroborated by a consideration of the red region in Figure 1 (recomputed for each of the 16 data), indicating likely models. The exponent pair of the Brody model was in the red region for all four gas car data. For the Gompertz model, it was close to the red regions of all six electric and hybrid car data and in addition of three gas car data (from AT and DE). For the West model (and the Bertalanffy model), it was in the red region of two conventional car data (gasoline-AT, diesel-DE) and three gas car data (from AT and DE). For logistic growth, it was in the red region of all electric and hybrid cars (except hybrid-DE).

Market Dynamics
We applied the BP models to explore the evolution of the market characteristics. We first obtained the best-fit BP models and then we used them to derive the LV Equation (2) for the market evolution. We selected gas cars as the outside goods as, for all countries, this was a rather marginalized drive technology. (For GE, we combined liquid and natural gas cars.) In Austria and Germany, the diesel market shares (computed from the best-fit growth model) decreased during the considered years 2011-2019 and 2012-2019. In Switzerland, diesel shares increased from 2005 to 2013 and decreased from 2013 to 2019. However, the loss of market shares did not necessarily mean a weaker position of diesel cars in terms of the market dynamics. This is because (disregarding the outside good), at the beginning, all technologies in all countries were in a state of competition, i.e., in Equation (2), all coefficients g i were positive. However, later, a predator-prey dynamic emerged, with diesel cars as the sole prey (g diesel < 0): this occurred in Germany and Switzerland in 2014 (t = 2.7 and t = 9.8, respectively) and in Austria in 2017 (t = 6.6). Another change occurred in 2018 (t = 13.2) in Switzerland, when gasoline cars became prey, too.

Did the Diesel Scandal Change the Market Dynamics?
We utilized model uncertainty to arrive at a more differentiated picture of market dynamics in terms of simulations, which augmented the outcomes for the best-fit models with information about other possible scenarios and their likelihood. As follows from our simulations (Table 3), in 2015 (year of the diesel scandal), the market position of diesel cars worsened in Germany and Switzerland. The role of diesel changed from competitor in 2013 and 2014 to (the sole) prey in 2015 and the following years. However, in Austria, other than as suggested by the best-fit models, the market standing of diesel most likely did not worsen, not even in the wake of the diesel scandal. Table 3 informs about the outcome of the simulations for the years 2013-2021 in more detail: for all countries, the scenario of pair-wise competition was most likely in 2013 and 2014, and the scenario of diesel as the sole prey was most likely for 2015-2019, except for Austria, where pair-wise competition remained the most likely scenario. Thus, for Germany and Switzerland, the simulations suggest that the best-fit model provided a premature date for the change in the role of diesel. However, in view of the high number of simulations with diesel as prey in 2013 (Germany) and 2014 (Germany and Switzerland), it was plausible that, in addition to the diesel scandal, another driver may have supported the change in 2015. Furthermore, other than for the best-fit model, the simulations suggested no change in the competitive pattern for Austria and no additional change after 2015 for Switzerland. Thus, also in this respect, the simulations corrected the best-fit model. This emphasizes the need for an analysis of model uncertainty to supplement the conclusions drawn from best-fit models.
We conclude that the worsening in Germany and Switzerland may not have been caused by the 2015 diesel scandal alone, as, for the simulations, already prior to the scandal, the standing of diesel cars deteriorated. We speculate that also government policies to support environmentally friendly cars by subsidies and privileges (e.g., free parking in cities) may have mattered. Such incentive-based policies [64] have been common for many years in several countries (e.g., Norway, The Netherlands, Sweden) and they were implemented also in Germany [65,66]. However, the VW scandal in turn may have influenced such policy responses, as, e.g., in Germany, it spurred a public debate about the prohibition of conventional cars, which has been recognized in other German-speaking countries as well [67]. Thus, the scandal may have eased the consensus finding in municipal councils for the introduction of disincentives directed specifically against diesel cars. For example, since 2018, the city of Stuttgart, Germany, has enforced a driving ban on old diesel vehicles to reduce air pollution [66]. Similar proposals were considered in Austria and Switzerland.

Forecasts
From the simulations for 2021, we note that the multi-model approach confirmed that modeling may not be the mathematical equivalent of the "crystal ball". For all countries, the attempt of a forecast for 2021 remained inconclusive. The simulations resulted in a proliferation of diverse scenarios and, at most, 50% of the simulations resulted in one of the three above-mentioned scenarios (which, in most previous years, accounted for 100%).
The reason is the broadening of the red bands for the extrapolation of the growth functions to the years 2020 and 2021 (as in Figures 2 and 4 for t = 9 and 10). Thus, for the present data, even forecasts for the near future that use likely models only carry significant model uncertainty. For this reason, in this paper, we do not consider forecasting. Furthermore, using trend analysis for this purpose assumes that the past trends will extend into the near future, while unexpected events (e.g., 2020 Covid-19 pandemic) may result in the falsification of such a prognosis.

Alternative Methodological Approaches
A simple alternative to our approach uses spline interpolation functions s i (t) of the market shares instead of the BP models. To define the LV model, it uses these splines in Equation (4) to define the coefficients of the LV Equation (2). While this approach reconstructs the market shares without errors, the splines result in complicated utility functions and consequently too many and rather questionable scenarios for the market dynamics (e.g., periods with all market goods in a state of symbiosis, the electric car as the sole predator, the hybrid car as the sole prey). By contrast, the trend analysis of this paper assumes more stability for the consumer preferences.
Our approach resulted in a surprisingly simple shape of the consumer utility functions (f i of Equation (4); c.f. Figure 5. This leads to the question of whether our modeling approach was overly complex. Why not start with the assumption that utility functions are parabolic in time? One study [68] has proposed this assumption, which leads to a simple alternative route to a LV model for the market shares, explained in Equation (9).
Thereby, z 0 is the remainder of the other z i to 100%. The starting point is the assumption that the utilities of the market goods are quadratic polynomials. The coefficients (a i , b i , c i ) of Equation (9) can be obtained by fitting the polynomial f i (t) to the observed utilities (they are computed from the observed market shares using formula (3)). Using Equation (9), the solutions z i of the LV Equation (2) are computed from the utilities [52,53]. Using this approach leads to similar results for the present data. The weaknesses of both approaches are similar, too: as has been illustrated [68] for different data, there may be large fluctuations in the observed market shares around the curves z i .
An advantage of the present method over the latter approach is the additional information (red band) that can be obtained from the multi-model approach. It allowed us to complement the conclusions from the best-fit model with simulations using models with a nearly best fit. We thereby could refine the market analysis from the best-fit model.

Discussion: Outlook on the Open Innovation in the Automotive Industry
The diesel scandal is a case study for the missed chances that arise if industry avoids open innovation. When VW and other companies realized that they could not fulfill the high environmental standards by conventional means (nitrogen oxide storage catalytic converters), they resorted to cheating, thereby exemplifying secrecy mentality at its worst. VW could have avoided this scandal if it had been more open to existing research on innovative technical solutions. For instance, new fuel injection technologies might resolve the pollution problems of diesel engines [69]. Fuel quality might be an issue, too [70]. Samsung is another example of a large corporation whose closed innovation strategy failed [71]: it had to withdraw its "Galaxy Note 7 smartphone" soon after the launch in 2016 (explosion of devices due to deficient batteries).
The resistance to open innovation in large corporations may be due to incompatible management structures [72]. However, there is pressure for a transformation of the traditional automotive industry (see below). As the example of Korean automakers shows, such a transformation of the automotive industry is already under way. They used open innovation in the form of "crowd innovation", which increased their performance, but at the cost of more patent disputes [73].
For a successful change, first, organizational, cultural, and sometimes also legal barriers need to be removed [74]. Amongst the cultural traits of firms that promote open innovation are a culture of failure acceptance, intra-venturing attitudes and intra-entrepreneurship, and a "transcendental-live firm culture" [75]. Such a culture is generally associated with companies from the information technology sector, e.g., Google or Apple. These are becoming new players in the automotive market, as software-related issues, such as autonomous driving or connected cars, are becoming ever more important [76,77]. Electric vehicles are contributing to this transformation, too. While the production of traditional cars requires high expertise in mechanics, electric cars can be produced by entrepreneurs, e.g., Tesla cars. Still, another pressure for the transformation of automotive industry comes from demands for a circular economy. For example, in 2015, the European Commission set up the action plan "Closing the Loop" [78]. Remanufacturing conventional cars from end-of-life parts saves energy, materials, and costs [79]. In developing countries, several companies specialize in this type of production.

Conclusions
We have studied the market diffusion of green cars in three countries. We started with the hypothesis that the diesel scandal might have led to increased sales of electric or hybrid cars. However, our paper has shown that, in the long run, the consumer behavior was surprisingly stable: the diesel scandal alone did not change the market characteristics. In Austria, diesel remained competitive and, in Germany and Switzerland, additional policy measures (prohibitions) were needed to change the consumer behavior, changing diesel from competitor to prey. However, as the car industry is currently undergoing a transformation towards more open innovation, we expect (from analogy to information technology) that the resulting changes from the supply side might create more demand for green cars. However, more research is needed to support this speculation.