Linear Regression Models for Time Series
The linear model
The simplest case is the simple linear regression model:
y_t = \beta_0 + \beta_1 x_t + \varepsilon_t
where:
- \beta_0 is the intercept: the predicted value of y when x = 0.
- \beta_1 is the slope: the average change in y for a one-unit change in x.
- \varepsilon_t is the error term: captures all other influences on y_t not explicitly modeled.
How are the \beta’s estimated? Ordinary Least Squares (OLS)
We cannot observe the \beta’s directly — we estimate them from data. The Ordinary Least Squares (OLS) method chooses the estimates \hat{\beta}_0, \hat{\beta}_1, \ldots, \hat{\beta}_k that minimize the sum of squared residuals:
\min_{\hat{\boldsymbol{\beta}}} \sum_{t=1}^{T} \hat{\varepsilon}_t^2 = \sum_{t=1}^{T} \left(y_t - \hat{\beta}_0 - \hat{\beta}_1 x_{1t} - \cdots - \hat{\beta}_k x_{kt}\right)^2
For the simple case (k = 1), taking partial derivatives and setting them to zero yields closed-form solutions:
\hat{\beta}_1 = \frac{\sum x_t y_t - T\bar{x}\bar{y}}{\sum x_t^2 - T\bar{x}^2}, \qquad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}
For the general case with k predictors, the solution extends naturally to matrix form:
\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}
where \mathbf{X} is the design matrix — a column of ones (for the intercept) plus one column per predictor — and \mathbf{y} is the vector of observed responses.
In practice, R handles all of this internally. But knowing the objective function matters: OLS penalizes large residuals more than small ones (because of the square), and under the classical assumptions the Gauss-Markov theorem guarantees the estimates are BLUE (Best Linear Unbiased Estimators) — regardless of whether we have one predictor or twenty.
Real examples
- Ice cream sales (y) and daily temperature (x).
- Nike revenue (y) and marketing spend (x).
- US consumption growth (y) and income growth (x).
- Mexico retail trade (y) and consumer confidence index (x).
Terminology: many names for y and x
The same concept appears under different names depending on the discipline:
| y (forecast variable) | x (predictor variables) |
|---|---|
| Dependent | Independent |
| Explained | Explanatory |
| Regressand | Regressor |
| Response | Stimulus / Covariate |
| Endogenous | Exogenous |
In time series forecasting, FPP3 uses forecast variable and predictor variables.
The origin of “regression” — Galton, 1886
The term was coined by Francis Galton while studying the relationship between parents’ height and their children’s height.
His finding: tall parents tend to have tall children, but on average their children are not as tall as them. Short parents tend to have children taller than themselves. There is a tendency to regress toward the mean.
Without this regression to the mean, the distribution of heights across generations would diverge — we would eventually have people of Hobbit stature and people of giant stature, with nothing in between.
Modern regression analysis generalizes this idea: studying the dependence of one variable on one or more others to predict its average value.
Are these models linear?
Exploratory analysis
Before fitting any model, we look at the data.
The scatter matrix shows correlations between all pairs of variables. Notice that some predictors are correlated with each other — we will return to this when we discuss multicollinearity.
Residual Diagnostics
Gauss-Markov assumptions and why they matter in time series
OLS produces BLUE estimators only if the classical assumptions hold. The one most frequently violated in time series is no autocorrelation in the errors:
\text{cov}(\varepsilon_t, \varepsilon_s) = 0 \quad \forall \; t \neq s
If residuals are autocorrelated, OLS estimates are still unbiased but no longer efficient — standard errors are wrong, p-values are wrong, and forecast intervals are wrong. This is why we always check the residuals.
- 1
- Extract the number of estimated parameters to use as degrees of freedom correction.
- 2
- A significant p-value indicates autocorrelated residuals — a sign the model is missing structure.
The residuals show clear autocorrelation. The simple model is not capturing all the dynamics of consumption. Let’s improve it.
Comparing simple vs. multiple
The multiple model substantially improves fit (\bar{R}^2 goes from ~0.15 to ~0.75) and the residuals look much closer to white noise.
The residuals vs. predictors plot is a useful diagnostic specific to multiple regression. If any panel shows a systematic pattern (curve, fan shape), it suggests a missing nonlinear term or interaction involving that predictor.
Three types of forecasts
Ex-ante forecasts use only information available at the time the forecast is made. The predictors must themselves be forecasted.
This is the “true” forecast — it is what you would actually produce in a real deployment.
# Step 1: forecast each predictor independently
us_change_fc_income <- us_change |>
model(ARIMA(Income)) |>
forecast(h = 8)
us_change_fc_prod <- us_change |>
model(ARIMA(Production)) |>
forecast(h = 8)
us_change_fc_sav <- us_change |>
model(ARIMA(Savings)) |>
forecast(h = 8)
us_change_fc_unemp <- us_change |>
model(ARIMA(Unemployment)) |>
forecast(h = 8)
# Step 2: construct future data with forecasted predictor values
us_change_future <- new_data(us_change, 8) |>
mutate(
Income = us_change_fc_income |> pull(.mean),
Production = us_change_fc_prod |> pull(.mean),
Savings = us_change_fc_sav |> pull(.mean),
Unemployment = us_change_fc_unemp |> pull(.mean)
)
# Step 3: forecast consumption using the projected predictors
us_change_fit_mult |>
forecast(new_data = us_change_future) |>
autoplot(us_change |> filter_index("2010 Q1" ~ .)) +
labs(title = "Ex-ante forecast: US Consumption",
y = "% change", x = NULL)
Ex-post forecasts use actual realized values of the predictors. The forecast variable y_t is still unknown, but we feed in the true future x_t values.
These are not “real” forecasts — they are used to evaluate the regression model in isolation, removing predictor forecast error from the evaluation.
# Use filter_index to create a train/test split
us_change_train <- us_change |>
filter_index(. ~ "2017 Q4")
us_change_test <- us_change |>
filter_index("2018 Q1" ~ .)
us_change_fit_expost <- us_change_train |>
model(
multiple = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
)
# Ex-post: we supply actual future predictor values
us_change_fit_expost |>
forecast(new_data = us_change_test) |>
autoplot(us_change |> filter_index("2014 Q1" ~ .)) +
labs(title = "Ex-post forecast: actual predictor values used",
y = "% change", x = NULL)
Scenario-based forecasts construct hypothetical future values for the predictors and ask: what would consumption look like under each scenario?
These are especially useful for stress-testing and strategic planning.
us_change_fit_scen <- us_change |>
model(TSLM(Consumption ~ Income + Savings + Unemployment))
us_change_scenarios <- scenarios(
optimistic = new_data(us_change, 4) |>
mutate(Income = c(0.5, 0.8, 0.6, 1.0),
Savings = c(0.1, -0.2, 0.1, -0.1),
Unemployment = -0.1),
pessimistic = new_data(us_change, 4) |>
mutate(Income = -0.5,
Savings = -0.4,
Unemployment = 0.2),
names_to = "Scenario"
)
us_change_fit_scen |>
forecast(new_data = us_change_scenarios) |>
autoplot(us_change |> filter_index("2014 Q1" ~ .)) +
labs(title = "Scenario-based forecast: US Consumption",
y = "% change", x = NULL)Beyond exogenous variables
TSLM() provides built-in terms that are particularly useful when working with time series:
trend()— a deterministic linear trend (t = 1, 2, \ldots, T).season()— seasonal dummy variables (one per season, with one omitted as baseline to avoid perfect collinearity).fourier(K)— Fourier terms for flexible seasonality (covered in Module 4).
We will illustrate with quarterly beer production in Australia.
Under the hood: seasonal dummies
Seasonal dummies are created automatically by season(), but it helps to see them explicitly:
Code
The dummy variable trap: if you create dummies for all m seasons and include an intercept, you have perfect multicollinearity. Always use m - 1 dummies. TSLM() handles this automatically.
There is a notable outlier in Q2 1994. We can handle it with an intervention variable.
Piecewise Linear Trends
When the trend changes direction or slope at one or more points, a piecewise linear model (also called a broken stick) can be appropriate.
We use the Boston Marathon winning times as an example — a series with clear structural breaks.
boston_fit <- boston_men |>
model(
linear = TSLM(Minutes ~ trend()),
exponential = TSLM(log(Minutes) ~ trend()),
piecewise = TSLM(Minutes ~ trend(knots = c(1940, 1980)))
)
boston_fc <- boston_fit |>
forecast(h = 10)
boston_men |>
autoplot(Minutes, color = "grey60") +
geom_line(data = fitted(boston_fit),
aes(y = .fitted, color = .model)) +
autolayer(boston_fc, alpha = 0.4, level = 95) +
labs(title = "Boston Marathon: three functional forms",
y = "Minutes", x = NULL, color = "Model") +
theme(legend.position = "top")- 1
- Log-transformed response: assumes exponentially declining times.
- 2
-
knotsspecify the breakpoints where the slope is allowed to change.
The piecewise model fits the historical data well, but notice what happens to the forecasts — particularly the linear and piecewise models. Always inspect the forecast, not just the in-sample fit. A model can fit the past perfectly and produce absurd predictions for the future.
