Data Science Fundamentals: Working with Time Series Data

Source: Google Gemini Nano Banana

Humans seemed to have an innate desire explain the future since the dawn of time. From early agriculture societies that used past experiences of the seasons to predict weather patterns, to attempts to read the stars through astrology, to the arcane methods of divination. But it took until the late 1800’s for this to develop into a formal practice, and even longer before the first time series models began to be applied. And if you’re someone who, like me, wants to dabble in the mystic arts of forecasting sans the woo-hoo elements, now is the best time to start honing your time series modeling skills. The era of big data has seen an explosion in the amount of information we have access to. And as the world becomes more chaotic, skills that can cut through the noise and achieve even a modicum of clarity become ever more valuable. So let’s get to learning!

Intro

For this exercise, we’ll be working in R and utilizing the tidyverts package for our time series applications. The main libraries we’ll be using with the tidyverts universe are fable and tsibble. This automatically loads tidyverse, so we’ll be utilizing that syntax as well.

The data we’ll be using is monthly oil and gas production in the US. The raw data comes from the Office of Natural Resources Revenue, a branch of the Department of the Interior, and can be found here. The available data dates back to 2003 and includes granular production details for all the US states. However, since our goal today is not to learn all there is to know about the oil and gas industry, nor how to clean a dataset, we’ll be using a curated version of this data from our good friends over at Kaggle. The final version contains close to 500k monthly observations dating from Jan 2015 to May 2025.

As is standard for machine learning, yₜ will refer to the observation at time t while ŷₜ refers to the forecast at time t.

Setup

The curated dataset contains a column for production disposition volumes with a cardinality over 30, so to simplify, we’ll aggregate the data to the monthly sum of oil production by State, and we’ll narrow our time series analysis to New Mexico as it is one of the largest oil producing states in the US.

oil_agg <- oil |>
  mutate(
    Date = mdy(Production.Date),
    Volume = parse_number(Volume)
  ) |> 
  filter(
    State != "Offshore",
    Disposition.Description %in% c(
      "Produced into Inventory Prior to Sales",
      "Sales-Royalty Due-MEASURED",
      "Sales-Royalty Due-Not Measured",
      "Sales-Royalty Not Due-MEASURED",
      "Sales-Royalty Not Due-FMP Not Assigned"
    )
  ) |> 
  select(-c(Land.Class, Land.Category, FIPS.Code, Offshore.Region, Production.Date, Disposition.Code, Disposition.Description)) |> 
  group_by(State, Commodity, Date) |> 
  summarise(
    Volume = sum(Volume, na.rm = TRUE),
    .groups = 'drop'
  )

Our next step is to convert this to what is known as a tsibble. In the tidyverts package, special versions of a data frame called a tsibble are used for time series modeling. They require a key, which are your time series objects, and an index as a unique identifier for each time series object.

oil_ts <- oil_agg |>
mutate(Date = yearmonth(Date)) |> 
  as_tsibble(
    key = c(State, Commodity),
    index = Date
  )

With this dataset, we have a total of 104 unique time series objects we could model (52 states and Oil and Gas commodities), indexed by month and year. The observations we are modeling are Volume (monthly production of the commodity).

Everything that follows you can replicate and do yourself by replicating the github repo. If you don’t want to deal with the code, then there is an app available to explore the data, build your own models and make forecasts. Links for both are below.

GitHub: Repo Link

App: Oil Exploration

Time Series Basics

When analyzing time series data, it’s helpful to think about the level (ℓₜ) of a time series. The level of a time series describes the overall baseline value of the series and is tied together by 4 main components.

• Trend: This component describes how the level changes over time. Do sales increase? Are click-through rates declining over time? What does the long-term path of the data look like?
• Cycle: These refer to irregular cycles that occur over long periods of time that you can’t set your clock to. Interest rate environments are a good example of this. Countries go through cycles of higher and lower interest rates, but the length and occurrence of these cycles is not fixed. This is often combined with the trend component to make up the trend-cycle.
• Season: This describes the regular fluctuations around the level that you can set your clock to. An obvious example is global temperatures following a regular seasonal pattern with peaks during the summer and troughs during the winter on repeat. A less obvious example is call-center volume. Call volume seasonal patterns can occur over hourly or weekly periods.
• Remainder: As the word suggests, this component contains anything not described by the other 3 components and often contains random, one-off events, like that day you got super motivated and worked out for 3 hours.

First Peek

With that in mind, let’s take a look at the time series plot of oil production in New Mexico.

Oil production in millions of barrels, where each barrel is 42 gallons. And no, it’s not measured in 42 gallons for that reason, it’s way more boring than that.

What do we see here? Is there a clear direction in oil production? Does there appear to be regular or irregular cycles in the data? What patterns can we pick out from this?

These questions are easier to answer by decomposing the data into it’s essential components. For time series data, an easy and robust way to do that is by modeling the data with Seasonal and Trend decomposition using Loess, or STL decomposition. This modeling can be done and visualized easily with tidyverts methods.

oil_model |> 
  filter(State == 'NM') |> 
  model(STL(`Oil (bbl)`, robust = F)) |> 
  components() |> 
  autoplot() +
  theme_minimal()

We first filter for the state of New Mexico then perform the modeling step. Robust is set to false meaning the algorithm uses standard least squares to fit the data making it more sensitive to outliers. The components method turns the results into a dable, which is tidyverts terminology for a decomposition data frame. Then we can see the results with autoplot which is a wrapper for ggplot.

These components change a lot by setting robust to TRUE, or by setting different lengths for the window argument in the trend and season parameters of the STL model. Play around with them yourself

The obvious component is the trend which is a steady line up and to the right. Less obvious from the time series plot is the seasonal component occurring on regular annual cycles. The variation in the seasons is also increasing with the overall level of the series. This is known as multiplicative behavior and will be important when we get to forecasting. The remainder component remains relatively steady until we get to mid-2020 to early 2022 when we see sharp swings in both directions (something crazy must have happened around then). The grey bars you see to the left of the plots are little helpers to get a sense of the scale of variation for each component. Each bar is the same length so the larger the bar, the smaller the overall variation of that component. Trend is the dominant component here with the remainder being a distant second followed up by season.

Modeling

Now that we know a little something about the data we’re working with, we can finally to get the fun part, forecasting. Last I checked, there was roughly a bajillion different models and algorithms to model and forecast time series data. Here we’ll compare two simple models, ETS and TSLM.

ETS stands for Error, Trend and Seasonality. This type of model is part of a broader family of Exponential Smoothing models which make forecasts based on weighted moving averages of previous observations. The more recent the observation, the higher the weight given. In other words, these models make the assumption that tomorrow is going to look more similar to today than last year.

In contrast to the statistical framework of the ETS models, TSLM, or Time Series Linear Models, assume that the observation yₜ can be expressed as a linear relationship with another time series xₜ. Being able to add external predictors to the model makes TSLM a bit trickier, but more interesting. However, at the end of the day, boring or not, results are what matter so let’s put the two head-to-head and see what happens.

The first step is to split our training data. For forecasting time series data, we can’t do a random shuffle as each observation is dependent on the previous one. We’ll keep it simple and use a cutoff date of 2024, meaning our training data will be from Jan 2015 to Dec 2023, and our validation data will be from Jan 2024 to May 2025.

The tidyverts framework makes it really easy to set up modeling, giving us more time to consider results and think of new ideas to try. We’ll use the below setup to compare the two models.

oil_fit <- oil_trn |> 
  filter(
    State == 'NM') |> 
  model(
    TSLM = TSLM(`Oil (bbl)` ~ Gas_Vol_Lag_one + Month + trend()),
    ETS = ETS(`Oil (bbl)`) 
  )

TSLM

So what’s going on here? For the TSLM model we see the variables Gas_Vol_Lag_one, Monthand trend(). Lets go over each.

• Gas_Vol_Lag_one: In life, it is often the case that events that happen today will affect what happens tomorrow, a month, or a year from now. So it is often helpful in time series modeling to include lagged versions of a variable that you think might have some predictive power on the target. In this case, we use the lag method to engineer a variable that lags gas production volume by one month so we can test if the previous months volume of gas production has any effect on this month's oil production.
• Month: There are two built in predictors worth adding to TSLM models and one of them is season(). In the tidyverts library, adding this argument to the TSLM model allows the algorithm to automatically detect the seasonal period of the time series, and use each period as a predictor. In monthly data like ours, that would show up in the model equation as season()year2 for February, seasonyear3() for March etc. Since we have the Month, we can use that to replace the season component and make the model equation more readable.
• trend(): This is the other built in parameter and it does what it says, it uses the trend component of the series to calculate forecasts.

ETS

The three components that make up an ETS model are error, trend and season. In the above set up, we leave it to the algorithm to automatically select the best components and their types. For a more thorough explanation of ETS models and each of their types, I’ll refer you to chapter 8 of Forecasting: Principles and Practice by Hyndman and Athanasopoulos. A simplified explanation is that the trend component can be additive or neutral, while error and season can be either of those or multiplicative. Multiplicative is appropriate when the component appears to be changing with the level of the series (told you that would be important later) and the coefficients of the components are expressed in relative terms instead of absolute like with additive.

After fitting the models, we can use report() and get a readout of the model coefficients and parameter values. I’ll put the important bits below, for the full model read out you can visit the app.

Model Reports

The ETS model chosen by the algorithm was ETS(M, Ad, M). This means the algorithm recognized the season changing in concert with the level of the series and the same for error. The little ‘d’ means it decided the trend component should be additive and damped. This provides an extra coefficient to the model that will flatten the trend line over time so that it doesn’t continue onward and upward for infinity. The values of the smoothing parameters are below.

alpha = 0.7423294 
beta  = 0.07237168 
gamma = 0.0001013897 
phi   = 0.9754635

Alpha (α) controls how the weights applied to previous observations change the further back you go. Higher levels give more weight to recent observations, indicating that oil production tomorrow will be more similar to how production is today than yesterday.

Beta (β) is the smoothing parameter for trend. Low values indicate a trend that isn’t changing often. In this case the model is recognizing the trend of oil production in New Mexico is fairly constant in one direction over time.

Gamma (γ) is the smoothing parameter for the season component. Very small values assume the seasonal component stays constant over time.

Phi (φ) controls how quickly the forecasted trend of production will flatten out. In the tidyverts library it’s constrained to being between 0.8 and 0.98.

For those familiar with the caret library, the TSLM model report looks very similar to the report output for logistic regression models. We get summary stats for the residuals, coefficient values, R2 and adjusted R2 and the f-stat. The interesting parts of this output are the p-values assigned to the predictors. In total, there are 14 predictors: the intercept, lagged gas production, trend, and 11 monthly indicator variables (one for each month, minus one, to avoid the dummy variable trap and maintain degrees of freedom). The only parameters shown as significant are the lagged gas production, trend, and the month of March. The adjusted R² is 0.9759. This is really high and is a good indicator of a superior model but it’s important to keep in mind that this measures how well the model fits historical data, not how well it predicts future data, which is what we’re hoping it will do.

Results

We can look at model reports all day, but models are only worth their salt if they do well on unseen data. We’ll look at a handful of metrics to get a handle on how well these models are doing.

oil_fit |> 
  forecast(new_data = oil_tst_nm) |> 
  accuracy(oil_tst_nm) |> 
  select(c(.model, ME, RMSE, MAE, MPE, MAPE))

We can use forecast to get predictions from our models. We can either use the argument h to state a certain period of time we want true forecasts for, or we can supply new_data to validate our models. Then accuracy gives us the performance of the models on several metrics and we can select the ones we care about most.

These look like big scary numbers but bear in the mind that New Mexico produces tens of millions of barrels each year. Root Mean Squared Error (RMSE) is a classic metric to compare models on because it penalizes large errors. The TSLM blows the ETS model out of the water on this one. Mean Error (ME) is always useful because it gives us an idea of the bias of the model. The high negative value for ETS indicates that it consistently over-forecasts the actual observation (ME in this instance is measured as yₜ - ŷₜ). The big numbers make it a bit challenging to tell what’s what, so Mean Percentage Error (MPE) and Mean Absolute Percentage Error (MAPE) can tell us the same things that ME and MAE do in relative terms. From the MAPE we can see that ETS is usually about 8% off and TSLM is about 4% off on average, twice as good as the ETS model.

They say a picture is worth a 1000 words, and that’s never more true than in data science.

Shaded regions represent prediction intervals (same thing as confidence intervals). Darker regions are an 80% confidence level, and lighter regions are 95%. The actual observations are represented by the black line

The most obvious observation is how much more uncertainty exists with the ETS model forecasts. Your instinct might be to treat that as a negative, but is it? In a perfect world we want to know the future with certainty, but in our extremely imperfect world, it’s better to predict the future realistically. The ETS model may be more realistic about the uncertainty of predicting oil production than the TSLM. A hint that that could be the case is the ETS model uncertainty expands the further out the forecast extends, unlike the TSLM uncertainty which stays the same. The future inherently becomes more unpredictable the further out you go, a simple fact reflected in the ETS model prediction intervals.

But this is the job of the data scientist no? The steps above are relatively simple, the hard part is making the decision that makes the most sense for your given situation. So what would you do? If the analysis stopped here, what decision would you make?

Me, I would drop the Month predictor and utilize the TSLM model. Other than March, none of the months showed up as significant to the model predictions, which lines up with the small seasonal component we saw in the decomposition. And whether or not the ETS model is more realistic about the uncertainty, the TSLM model did way better on the validation data.

Thankfully for us, the analysis doesn’t have to stop here! These are preliminary steps, and there are other diagnostics we can use to assess our models, such as autocorrelation plots, other techniques to assess performance over multiple time periods like cross-validation, and other models to try like ARIMA, all of which we’ll get into next time. Until then, stay curious!

References

Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. Forecasting: Principles and Practice (3rd ed)

O’Hara-Wild, M., Hyndman, R., & Wang, E. (2024). fable: Forecasting Models for Tidy Time Series (R package version 0.4.1). https://fable.tidyverts.org

Topuz, P. (2025). US Oil & Gas Production & Disposition 2015–2025 [Data set]. Kaggle. https://www.kaggle.com/datasets/pinuto/us-oil-and-gas-production-and-disposition-20152025

Wang, E., & Hyndman, R. J. (2021). tsibble: Tidy Data Structures for Time Series (R package version 1.1.6). https://tsibble.tidyverts.org