No change

The simplest model is to assume no change. Hot or cold, we assume to sale the same amount of ice cream. I can express this as an ordinary differential equation (ODE):

\[ \frac{d\mu(t)}{dt} = 0,\quad\mu(0) = \mu_0 \] It’s solution is \(\mu(t) = \mu_0\) and hence, we could assume for the observable data \(y(t)\): \[ y(t) \sim \mathcal{N}(\mu_0, \sigma^2) \]

Of course, this can’t be right. I would expect higher sales of ice cream for higher temperatures.

Constant change

By assuming constant change I mean that the increase in sales is constant, regardless of the starting temperature, i.e. a change in temperature from \(7^\circ C\) to \(12^\circ C\) will cause the same increase in sales as from \(22^\circ C\) to \(27^\circ C\), or to put this into a ODE:

\[ \frac{d\mu(t)}{dt} = b,\quad \mu(0)=a \] Solving the ODE gives us the classic linear equation for the mean: \[ \mu(t) = a + b\,t \] Now, as we are assuming constant change, I should also assume that the variance \((\sigma^2)\) around my data doesn’t change either and hence, I could model my data as Gaussian again:

\[ y(t) \sim \mathcal{N}(\mu(t), \sigma^2) = \mathcal{N}(a + b\,t, \sigma^2) \]

This model makes a little more sense, but it could also predict negative expected sales, if \(a\) is negative. In addition, generating data from this distribution will produce negative values occasional as well, even if \(a>0\).

Constant sales growth rate

Assuming a constant sales growth rate means that a \(1^\circ C\) increase in temperature is responded by the same proportional change in sales, i.e. for every \(1^\circ C\) increase the sales would increase by \(b\)%. Thus, we are moving from an additive to a multiplicative scale. We can describe the change in the median sales as: \[ \frac{d\mu(t)}{dt} = \mu(t) \, b,\quad \mu(0)=a \] The solution to this ODE is exponential growth (\(b>0\)), or decay (\(b<0\)): \[ \mu(t) = a \, \exp(b\,t) \] On a log-scale we have: \[ \log{\mu(t)} = \log{a} + b\,t \]

This is a linear model again. Thus, we can assume constant change and variance \((\sigma^2)\) on a log-scale. A natural choice for the data generating distribution would be:

\[ \log{y}(t) \sim \mathcal{N}(\log{a} + b\,t, \sigma^2) \]

Or, in other words, I believe the data follows a log-normal distribution. The expected sale for a given temperature is \(E[y(t)]=\exp{(\mu(t)+\sigma^2/2)}\), its median is \(\exp{(\mu(t))}\). In addition this also means that we assume a constant coefficient of variation \((CV)\) on the original scale:

\[ \begin{aligned} CV &=\frac{\sqrt{Var[y]}}{E[y]} \\ & = \frac{\sqrt{(e^{\sigma^2} - 1)e^{2\mu(t)+\sigma^2}}}{e^{\mu(t)+\sigma^2/2}} \\ & = \sqrt{e^{\sigma^2} - 1} \\ & = const \end{aligned} \]

Although this model avoids negative sales, it does assume exponential growth as temperatures increase. Perhaps, this a little too optimistic.

Constant demand elasticity

A constant demand elasticity of \(b = 60\)% would mean in this context that a 10% increase in the temperature would lead to a 6% increase in demand/ sales.

Thus, the relative change in sales is proportional to the relative change in temperature:

\[ {\frac{d\mu(t)}{\mu(t)}} = b {\frac{dt}{t}} \] Re-arranging and adding an initial value gives us the following ODE:

\[ \frac{d\mu(t)}{dt} = \frac{b}{t}\mu(t),\quad\mu(1)=a \]

It’s solution is:

\[ \mu(t) = \exp{(a + b \log(t))} = \exp{(a)} \cdot t^b \]

Now take the log on both sides and we are back to linear regression equation:

\[\log(\mu(t)) = a + b \log(t)\] Thus, we are assuming constant change on a log-log scale. We can assume a log-normal data generating distribution again, but this time we measure temperature on a log-scale as well.

\[ \log{y(t)} \sim \mathcal{N}(a + b\log{t}, \sigma^2) \]

This model, as well as the log-transformed model above, predicts ever increasing ice cream sales as temperature rises.

Constant market size

However, if I believe that my market/ sales opportunities are limited, i.e. we reach market saturation, then a logistic growth model would make more sense.

Let’s assume as temperature rises, sales will increase exponential initially and then tail off to a saturation level \(M\). The following ODE would describe this behaviour:

\[ \frac{d\mu(t)}{dt} = b\, \mu(t) (1 - \mu(t)/M),\quad\mu(0)=\mu_0 \]

Integrating the ODE provides us with the logistic growth curve (often used to describe population growth):

\[ \mu(t) = \frac{ M }{1 + \frac{M-\mu_0}{\mu_0}e^{-b\,t}} = \frac{\mu_0 M e^{b\,t}}{\mu_0 (e^{b\,t} - 1) + M} \]

Setting \(M=1\) gives me the proportion of maximum sales, and setting \(a:=(1-\mu_0)/\mu_0\) gives us:

\[ \begin{aligned} \mu_p(t) &= \frac{ 1 }{1 + \frac{1-\mu_0}{\mu_0}e^{-b\,t}} \\ & = \frac{1}{1 + e^{-(a+b\,t)}} \end{aligned} \]

Nice, this is logistic regression. Defining \(\mbox{logitic}(u) := \frac{1}{1 + e^{-u}}\), we can see that the above has a linear expression in \(u\).

Thus, if I set \(z(t):=y(t)/M\), I could assume the Binomial distribution as the data generating process:

\[ z(t) \sim \mbox{Bin}\left(M, \mbox{logitic}(\mu_p(t))\right) \]

But I don’t know \(M\). Thus, I have to provide a prior distribution, for example a Poisson distribution with parameter (my initial guess) \(M_0\)

\[ M \sim \mbox{Pois}(M_0) \]

Or, alternatively I could consider the Poisson distribution as the data generating distribution, with its mean and variance being described by the non-linear growth curve above (this is not a GLM, unlike the models so far):

\[ y(t) \sim \mbox{Pois}(\mu(t)) = \mbox{Pois}\left(\frac{\mu_0 M e^{b\,t}}{\mu_0 (e^{b\,t} - 1) + M}\right) \]

The Poisson distribution has the property that the mean and variance are the same. If \(y \sim \mbox{Pois}(\lambda)\), then \(E[y] =\mbox{Var}[y] = \lambda\), but more importantly, we are generating integers, in line with the observations.

So, let’s generate some data to check if this approach makes sense in the temperature range from \(-5^\circ C\) to \(35^\circ C\).

Let’s put priors over the parameters: \[ \begin{aligned} \mu_0 &\sim \mathcal{N}(10, 1),\quad \mbox{ice creams sold at }0^\circ C\\ b & \sim \mathcal{N}(0.2, 0.02^2),\quad \mbox{inital growth rate} \\ M & \sim \mathcal{N}(800, 40^2) ,\quad \mbox{total market size} \end{aligned} \]

temperature <- c(-5:35)
n <- length(temperature)
nSim <- 10000
set.seed(1234)
# Priors
mu0 = rep(rnorm(nSim, 10, 1), each=n)
b = rep(rnorm(nSim, 0.2, 0.02), each=n)
M = rep(rnorm(nSim, 800, 40), each=n)
lambda = matrix(M/(1 + (M - mu0)/mu0 * exp(-b * rep(temperature, nSim))),
                ncol=n)

With those preparations done I can simulate data:

y_hat_p <- matrix(rpois(nSim*n, lambda), nrow=n)
probs <- c(0.025, 0.25, 0.5, 0.75, 0.975)
pp=t(apply(y_hat_p, 1, quantile, probs=probs))
matplot(temperature, pp, t="l", bty="n", 
        main=paste("Poisson percentiles:",
                   "2.5%, 25%, 50%, 75%, 97.5%", sep="\n"))

This looks reasonable, but the interquartile range looks quite narrow.

Hence, I want to consider the Negative-Binomial distribution, which has an additional shape parameter \(\theta\) to allow for over-dispersion, i.e. the variance can be greater than the mean. A Negative Binomial distribution with parameters \(\mu\) and \(\theta\) has mean \(\mu\) and variance \(\mu + \mu^2/\theta\). With \(\theta \to \infty\) we are back to a Poisson distribution.

library(MASS) # provides negative binomial
y_hat_nb <- matrix(rnegbin(nSim*n, lambda, theta=20), nrow=n)
nbp <- t(apply(y_hat_nb, 1, quantile, probs=probs))
matplot(temperature, nbp, t="l", bty="n", 
        main=paste("Negative Binomial percentiles:",
        "2.5%, 25%, 50%, 75%, 97.5%", sep="\n"))

This looks a little more realistic and adds more flexibility to my model.

Finally, I am getting to a stage where it makes sense to make use of my data.