Modeling a Short Time Series
Imagine the following scenario. Your company has launched a new division whose task it is to sell bikes in New York City. After observing sales for the first quarter since launch, it is time to forecast sales for the holidays and the upcoming year. The following graph shows the daily sales so far.
The data we use in the following is publicly available Citi Bike data from station 360 in New York City. You can access the raw data here.
A Naive Benchmark Model
Give someone such a time series and they tend to grab ready-made solutions such as Hyndman’s forecast or Sean J. Taylor’s and Ben Letham’s prophet package (for good reason). Let’s see how far these methods get us.
sales <- hciti %>%
pull(rides) %>%
ts(frequency = 7)
arima_fit <- sales %>%
forecast::auto.arima()
arima_fc <- forecast::forecast(arima_fit, h = 180)
While it is super simple to fit an ARIMA model, it’s often not sufficient to throw the automatic procedure at the data. Here, we quickly lose the weekly seasonality of the data in our forecasts.
The first trial can be improved with a little bit of intervention.
xreg_fourier <- forecast::fourier(sales, K = 3)
xreg_fourier_future <- forecast::fourier(sales, K = 3, h = 180)
arima_fit <- forecast::auto.arima(sales, seasonal = FALSE,
xreg = xreg_fourier)
arima_fc <- forecast::forecast(arima_fit, h = 180,
xreg = xreg_fourier_future)
While this second iteration might be a better model of the weekly seasonality in our data, it fails to incorporate most of our problem-specific information.
Let’s take a step back and actually think about what we’re trying to do.
Existing Business Knowledge
First of all, we are dealing with count data: We sell 1 bike, we sell 102 bikes, but we don’t sell neither 8.1 bikes nor -15 bikes. As we can see in the plots above, this knowledge is not well represented by the ARIMA models with its assumption of Normal-distributed errors. The prediction intervals quickly move into the negative area which is not sensible and implies that we attribute not enough probability to the positive outcomes.1 One can try to account for this via a Box-Cox transformation, such as forecast::auto.arima(sales, seasonal = FALSE, xreg = xreg_fourier, lambda = 0).
Second, the forecast does not incorporate a yearly seasonality. Even though we have not observed a year of sales, our boss might have the strong belief that bike sales have a seasonality following nature’s seasons, just like the demand for ice cream: Higher in summer, lower in winter. The problem with any kind of seasonality, though, is that common advice is to only model the seasonality once you have at least two periods of data. That is less of a problem with the weekly seasonality, as we already observed several weeks. But waiting for two years before we can even try to incorporate yearly seasonality? That’s going to leave our stakeholders unsatisfied.
This is what happens when you model yearly seasonality naively with a quarter of data using maximum likelihood.
The most frustrating part about not being able to model the long seasonality is that we as data scientists as well as our colleagues know that crucial knowledge about the business is being neglected. How can we expect anyone to trust the forecast when it ignores obvious business knowledge?
Another expectation from the business could be to see growth in the demand for our bikes over time. Even if there is no growth to be seen yet, we might have to account for growth because of planned marketing, because the first spring and summer are coming up, or because we previously observed growth with other products.
Hell breaks loose when we add another free parameter for a linear trend to the yearly seasonality.
The results from Facebook’s prophet package expose similar difficulties.
Collecting Everything for a Data Generating Process
In the previous paragraphs, we loosely specified assumptions we and our colleagues would like to see incorporated in a forecast of the bike sales. We also saw that we cannot apply a typical time series model on the few data we have and expect a good result that is in line with our knowledge. Thus, we will now try to build a probabilistic model that let’s us include the prior knowledge. Let us formalize the assumptions stated so far.
We model random variables \(Y_t\) representing sales with a count distribution such as Poisson or a Negative Binomial distribution. For the rest of the discussion, we continue with a Negative Binomial parameterized by \(\mu_t\) and \(\phi\), such that the mean and variance of \(Y_t\) are given by \(\text{E}[Y_t] = \lambda_t\) and \(\text{Var}[Y_t] = \lambda_t + \frac{\lambda_t^2}{\phi}\):
\[ Y \sim \text{NegBin}(\lambda_t, \phi), \] where \(\lambda_t, \phi > 0\).
Since we’re dealing with time series, we model the mean with a damped dynamic that combines information from the previous observation \(y_{t-1}\), the latent mean at the previous time step \(\lambda_{t-1}\), as well as the general level \(\mu_t\):
\[ \lambda_t = (1 - \delta - \gamma ) \cdot \mu_t + \delta \cdot \lambda_{t-1} + \gamma \cdot y_{t-1}, \] where \(0 \leq \delta \leq 1\) and \(0 \leq \gamma \leq 1-\delta\).
We model the general level \(\mu_t\) of the time series as a regression using dummy variables \(D_i\) for every day, as well as Fourier terms for the long yearly seasonality.2 See this blog post by Rob J. Hyndman or this paper on Fourier Recurrent Units in LSTMs for more information and motivation. Given that the mean of a Negative Binomial has to be positive, we use a log-link function:
Fourier series terms of different order for a seasonality of period 365.25.
\[ \log(\mu_t) = \alpha_0 + \sum_{i = 1}^6 \alpha_i D_i + \tau \frac{t}{m} + \sum_{k=1}^K \Big(\beta_k \sin \big(\frac{2\pi kt}{m} \big) + \gamma_k \cos \big(\frac{2\pi kt}{m} \big)\Big) \]
where \(m = 365.25\). \(K\) defines how complex we let the Fourier terms and corresponding seasonal pattern be; given that we are mostly interested in a cycle corresponding to the natural seasons, \(K = 6\) should be more than sufficient. The \(\tau \frac{t}{m}\) term models a linear trend of the log-transformed mean.
To make the data generating process complete, we add prior distributions to model the uncertainty for every parameter.
Following the Prior Choice Recommendations for the probabilistic programming language Stan, we don’t put a prior on \(\phi\) directly, but instead on the transformed \(\hat{\phi} = \frac{1}{\sqrt{\phi}}\) with
\[ \hat{\phi} \sim N(0,0.5). \]
As priors on the lagged mean and observation components, we use Gamma priors to enforce positive values with \(\delta \sim \text{Gamma}(1,10)\) and \(\gamma \sim \text{Gamma}(0.5,10)\).
For the coefficients in the regression equation, we don’t need any transformation. Given that we have quite some data for every day and since they are not the focus of this article, we give weakly-informative priors for every day.
\[ \alpha_i \sim N(0,0.5). \]
The prior for the intercept corresponds to the level of sales we roughly expect. Note that we model the intercept on the logarithmic scale: We expect around \(\exp(4) \approx 55\) sales per day.
\[ \alpha_0 \sim N(4,0.5) \]
More interesting, and the actual focus of this article, however, are the priors on the trend and Fourier terms coefficients. We could, of course, choose penalizing priors to account for the fact that we have not observed much data. We might not want to estimate long-term growth based on the first few weeks of data. Similar for the Fourier terms; we have seen above in the ARIMA examples that the parameters add too much flexibility as that they could be estimated reliably from our amount of data.
Instead, we want to very much inform the model about what parameters it can expect in the real world. Having weakly informing priors does not make sense; we have a clear idea of how much growth makes sense, and what the seasonality might look like. Let’s start with the trend.
Choosing a Prior for the Trend
One advantage of the log-link parameterization of the mean is that we may interpret the coefficients as follows: A one unit change in a given variable corresponds to a change in \(\mu_t\) corresponding to 100% multiplied by the variable’s coefficient, holding every other variable constant. So for the trend as defined above, \(\tau\) gives the percentage change in \(\mu_t\) from one year to the next. Thus, if the business very much plans with a growth of, let’s say, 3% year-over-year, then we should choose \(\text{E}[\tau] = 0.03\).
If we additionally believe very much in slight positive growth, then a prior such as
\[ \alpha_0 \sim N(0.03, 0.02) \]
might be reasonable. Note that this does not imply negative growth is impossible—it just means we think it’s unlikely. Indeed, our prior quantifies a negative growth scenario with a probability of about 6.68%.
Choosing a Prior for the Yearly Seasonality
Thinking about a useful linear combination of sine and cosine curves that combined portray the right amount of up and down over a year is a little more involved than choosing a single YoY growth. It would surely be possible to come up with reasonable coefficients for the Fourier terms by thinking hard or playing around, but we don’t really have an intuition on what reasonable priors on the coefficients look like. We do, however, have an idea about the functions that these priors induce together. So we will try a different approach here.
Three years of daily maximum temperatures for New York City. 
Posterior linear combination of the Fourier terms for the temperature model. 
Posterior mean linear combination of the Fourier terms for the temperature model. 