library(TSstudio)
data(USgas)
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")40 Forecasting with Linear Regression
The linear regression model is one of the most common methods for identifying and quantifying the relationship between a dependent variable and a single (univariate linear regression) or multiple (multivariate linear regression)independent variables. This model has a wide range of applications, from causal inference to predictive analysis and, in particular, time series forecasting.
40.1 The Linear Regression
The primary usage of the linear regression model is to quantify the relationship between the dependent variable \(Y\) (also known as the response variable) and the independent variable/s \(X\) (also known as the predictor, driver, or regressor variables) in a linear manner. In other words, the model expresses the dependent variable as a linear combination of the independent variables.
A linear relationship between the dependent and independent variables can be generalized by the following equations:
- In the case of a single independent variable, the equation is as follows:
Here we are modeling three of the four time series components: level (\(\beta_0\)), trend (\(\beta_1\)) and noise (\(\epsilon\)). Seasonality is not modeled.
- For n independent variables, the equation looks as follows:
The term linear, in the context of regression, referred to the model coefficients, which must follow a linear structure (as this allows us to construct a linear combination from the independent variables). On the other hand, the independent variables can follow both a linear and nonlinear formation.
A linear regression model provides an estimation for those coefficients (that is, \(\beta_0\), \(\beta_1\),…, \(\beta_n\)), which can be formalized by the following equations:
- For the univariate linear regression model, the equation is as follows:
- For the multivariate linear regression model, the equation is as follows:
40.1.1 Estimation of coefficients
The estimation of the model’s coefficients is based on the following two steps:
Define a cost function (also known as loss function) —setting some error metric to minimize
Apply mathematical optimization for minimizing the cost function.
The most common estimation approach is applying the Ordinary Least Squares (OLS) method as an optimization technique for identifying the coefficients that minimize the residuals sum of squares:
\[\sum_{i=1}^{N} \hat{\epsilon_i^2}\]
Besides the \(OLS\), there are other multiple estimation methods, such as the maximum likelihood, method of moments, and Bayesian.
40.2 Forecasting with Linear Regression
The linear regression model, unlike the traditional time series models such as the \(ARIMA\) or \(Holt-Winters\), was not designed explicitly to handle and forecast time series data. Instead, it is a generic model with a wide range of applications from causal inference to predictive analysis. Therefore, forecasting with a linear regression model is mainly based on the following two steps:
Identifying the series structure, key characteristics, patterns, outliers, and other features.
Transforming those features into input variables and regressing them with the series to create a forecasting model.
40.2.1 Forecasting the trend and seasonal components
For the sake of simplicity, we can drop the cycle component as it typically merged into the trend component (or ignored the cycle component). We can now transform the (additive and multiplicative) time series decomposition equations for a linear regression model:
Where:
\(Y_t\) represents a time series with \(n\) observations.
\(T\), an independent variable with \(n\) observations, represents the series trend component.
\(S_t\), an independent variable with \(n\) observations, represents the series seasonal component.
\(\epsilon_t\), the regression error term, represents the irregular component or any pattern that is not captured by the series trend and seasonal component.
\(\beta_0, \beta_1, \beta_2\) represent the level (model intercept), and coefficients of the trend and seasonal components, respectively.
For the sake of convenience and in the context of working with time series, we will change the observations notation from \(i\) to \(t\), as it represents the time dimension of the series.
The transformation of a series with a multiplicative structure into a linear regression is done in two steps:
- The transformation of the series by applying a
logtransformation for both sides of the equations of the series:
- Once the series is transformed into an additive structure, the transformation to a linear regression formation:
40.2.2 Features engineering of the series components
Before creating the regression inputs that represent the series trend and seasonal components, we first need to understand their structure.
Example:
As you can see in the series plot, and as we saw on the previous chapters, USgas is a monthly series with a strong monthly seasonal component and fairly stable trend line.
ts_info(USgas) The USgas series is a ts object with 1 variable and 238 observations
Frequency: 12
Start time: 2000 1
End time: 2019 10 We can explore the series components structure with the ts_decompose function further:
ts_decompose(USgas)You can see in the preceding plot that the trend of the series is fairly flat between 2000 and 2010, and has a fairly linear growth moving forward. Therefore, the overall trend between 2000 and 2018 is not strictly linear. This is an important insight that will help us to define the trend input for the regression model.
Next we have to transform the series from a ts object to a data.frame object. Therefore, we will utilize the ts_to_prophet function from the TSstudio package:
USgas_df <- ts_to_prophet(USgas)
head(USgas_df) ds y
1 2000-01-01 2510.5
2 2000-02-01 2330.7
3 2000-03-01 2050.6
4 2000-04-01 1783.3
5 2000-05-01 1632.9
6 2000-06-01 1513.1Now we can start to create the regression input features. The first feature we will create is the series trend. A basic approach for constructing the trend variable is by indexing the series observations in chronological order:
USgas_df$trend <- 1:nrow(USgas_df)The second feature we want to create is the seasonal component. Since we want to measure the contribution of each frequency unit to the oscillation of the series, we will use a categorical variable for each frequency unit. Therefore, we will create a categorical variable with 12 categories, each category corresponding to a specific month of the year in the USgas series.
library(lubridate)
USgas_df$seasonal <- factor(month(USgas_df$ds, label = T), ordered = FALSE)
head(USgas_df) ds y trend seasonal
1 2000-01-01 2510.5 1 ene
2 2000-02-01 2330.7 2 feb
3 2000-03-01 2050.6 3 mar
4 2000-04-01 1783.3 4 abr
5 2000-05-01 1632.9 5 may
6 2000-06-01 1513.1 6 junLast but not least, before we start to regress the series with those features, we will split the series into a training and testing partition. We will set the last 12 months of the series as a testing partition:
h <- 12 # setting a testing partition length
train <- USgas_df[1:(nrow(USgas_df) - h), ]
test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]Now we can proceed and review how the regression model captures each one of the components separately and all together.
40.2.3 Modeling the series trend and seasonal components
TREND COMPONENT
We will first model the series trend by regressing the series with the trend variable, on the training partition:
Example:
md_trend <- lm(y ~ trend, data = train)
summary(md_trend)
Call:
lm(formula = y ~ trend, data = train)
Residuals:
Min 1Q Median 3Q Max
-547.2 -307.4 -153.2 333.1 1052.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1751.0074 52.6435 33.26 < 2e-16 ***
trend 2.4489 0.4021 6.09 4.86e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 394.4 on 224 degrees of freedom
Multiple R-squared: 0.1421, Adjusted R-squared: 0.1382
F-statistic: 37.09 on 1 and 224 DF, p-value: 4.861e-09As you can see from the preceding regression output, the coefficient of the trend variable is statistically significant to a level of 0.001. However, the adjusted R-squared of the regression is fairly low, which generally makes sense, as most of the series variation of the series is related to the seasonal pattern as we saw in the plots previously.
Explanation: As you can note from the preceding regression output, the fourth column represents the p-value of each one of the model coefficients. The p-value provides the probability that we will reject the null hypothesis given it is actually true, or the type I error. Therefore, for the p-value smaller than \(\alpha\), the threshold value, we will reject the null hypothesis with a level of significance of \(\alpha\), where typical values of \(\alpha\) are 0.1, 0.05, 0.01, and so on.
Our next step is to use the model we created to predict the fitted values on the training partition and the forecasted values on the testing partition.
train$yhat <- predict(md_trend, newdata = train)
test$yhat <- predict(md_trend, newdata = test)We will create a utility function that plots the series and the model output, utilizing the plotly package:
library(plotly)
plot_lm <- function(data, train, test, title = NULL){
p <- plot_ly(data = data,
x = ~ ds,
y = ~ y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = ~ train$ds,
y = ~ train$yhat,
line = list(color = "red"),
name = "Fitted") %>%
add_lines(x = ~ test$ds,
y = ~ test$yhat,
line = list(color = "green", dash = "dot", width = 3),
name = "Forecasted") %>%
layout(title = title,
xaxis = list(title = "Year"),
yaxis = list(title = "Billion Cubic Feet"),
legend = list(x = 0.05, y = 0.95))
return(p)
}Let’s set the inputs of the plot_lm function with the model output:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend Component of the Series")Overall, the model was able to capture the general movement of the trend, yet a linear trend may fail to capture the structural break of the trend that occurred around 2010.
Last but not least, for comparison analysis, we want to measure the model error rate both in the training and the testing sets:
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_trend[1] 0.1644088 0.1299951SEASONAL COMPONENT
The process of modeling and forecasting the seasonal component follows the same process as we applied with the trend, by regressing the series with the seasonal variable we created before:
Example:
md_seasonal <- lm(y ~ seasonal, data = train)
summary(md_seasonal)
Call:
lm(formula = y ~ seasonal, data = train)
Residuals:
Min 1Q Median 3Q Max
-608.98 -162.34 -50.77 148.40 566.89
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2774.28 49.75 55.759 < 2e-16 ***
seasonalfeb -297.92 70.36 -4.234 3.41e-05 ***
seasonalmar -479.10 70.36 -6.809 9.77e-11 ***
seasonalabr -905.28 70.36 -12.866 < 2e-16 ***
seasonalmay -1088.42 70.36 -15.468 < 2e-16 ***
seasonaljun -1105.49 70.36 -15.711 < 2e-16 ***
seasonaljul -939.35 70.36 -13.350 < 2e-16 ***
seasonalago -914.12 70.36 -12.991 < 2e-16 ***
seasonalsep -1114.74 70.36 -15.843 < 2e-16 ***
seasonaloct -1022.21 70.36 -14.527 < 2e-16 ***
seasonalnov -797.53 71.33 -11.180 < 2e-16 ***
seasonaldic -256.67 71.33 -3.598 0.000398 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 216.9 on 214 degrees of freedom
Multiple R-squared: 0.7521, Adjusted R-squared: 0.7394
F-statistic: 59.04 on 11 and 214 DF, p-value: < 2.2e-16Since we regress the dependent variable with a categorical variable, the regression model creates coefficients for 11 out of the 12 categories, which are those embedded with the slope values. As you can see in the regression summary of the seasonal model, all the model’s coefficients are statistically significant. Also, you can notice that the adjusted R-squared of the seasonal model is somewhat higher with respect to the trend model (0.78 as opposed to 0.1).
Before we plot the fitted model and forecast values with the plot_lm function, we will update the values of yhat with the predict function:
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)Now we can use the plot_lm function to visualize the fitted model and forecast values:
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Seasonal Component of the Series")As you can see in the preceding plot, the model is doing a fairly good job of capturing the structure of the series seasonal pattern. However, you can observe that the series trend is missing. Before we add both the trend and the seasonal components, we will score the model performance:
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_seasonal[1] 0.08628973 0.21502100The high error rate on the testing set is related to the trend component that was not included in the model.
The next step is to join the two components into one model and to forecast the feature values of the series:
md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)
Call:
lm(formula = y ~ seasonal + trend, data = train)
Residuals:
Min 1Q Median 3Q Max
-514.73 -77.17 -17.70 85.80 336.95
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2488.8994 32.6011 76.344 < 2e-16 ***
seasonalfeb -300.5392 41.4864 -7.244 7.84e-12 ***
seasonalmar -484.3363 41.4870 -11.674 < 2e-16 ***
seasonalabr -913.1334 41.4880 -22.010 < 2e-16 ***
seasonalmay -1098.8884 41.4895 -26.486 < 2e-16 ***
seasonaljun -1118.5855 41.4913 -26.960 < 2e-16 ***
seasonaljul -955.0563 41.4936 -23.017 < 2e-16 ***
seasonalago -932.4482 41.4962 -22.471 < 2e-16 ***
seasonalsep -1135.6874 41.4993 -27.366 < 2e-16 ***
seasonaloct -1045.7687 41.5028 -25.198 < 2e-16 ***
seasonalnov -808.0016 42.0617 -19.210 < 2e-16 ***
seasonaldic -269.7642 42.0635 -6.413 9.05e-10 ***
trend 2.6182 0.1305 20.065 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 127.9 on 213 degrees of freedom
Multiple R-squared: 0.9142, Adjusted R-squared: 0.9094
F-statistic: 189.2 on 12 and 213 DF, p-value: < 2.2e-16Regressing the series with both the trend and the seasonal components together provides additional lift to the adjusted R-squared of the model from 0.78 to 0.91. This can be seen in the plot of the model output:
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend and Seasonal Components of the Series")Let’s measure the model’s \(MAPE\) score on both the training and testing pPartitions:
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1[1] 0.04774945 0.09143290Regressing the series with both the trend and the seasonal components provides a significant lift in both the quality of fit of the model and with the accuracy of the model.
CAPTURING A NON-LINEAR TREND
When looking at the plot of the model fit and forecast, you can notice that the model trend is too linear and missing the structural break of the series trend. This is the point where adding a polynomial component for the model could potentially provide additional improvement for the model accuracy.
A simple technique to capture a non-linear trend is to add a polynomial component to the series trend in order to capture the trend curvature over time. We will use the I argument, which allows us to apply mathematical operations on any of the input objects. Therefore, we will use this argument to add a second degree of the polynomial for the trend input:
Example:
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)
summary(md2)
Call:
lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
Residuals:
Min 1Q Median 3Q Max
-468.47 -54.66 -2.21 63.11 294.32
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
seasonalfeb -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
seasonalmar -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
seasonalabr -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
seasonalmay -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
seasonaljun -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
seasonaljul -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
seasonalago -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
seasonalsep -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
seasonaloct -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
seasonalnov -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
seasonaldic -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
trend -1.270e+00 4.472e-01 -2.840 0.00494 **
I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 109.1 on 212 degrees of freedom
Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16Adding the second-degree polynomial to the regression model did not lead to a significant improvement of the goodness of fit of the model.
On the other model, as you can see in the following model output plot, this simple change in the model structure allows us to capture the structural break of the trend over time:
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)
plot_lm(data = USgas_df,
train = train,
test = test,
title = "Predicting the Trend (Polynomial) and Seasonal Components of the Series")As we can see from the model following the MAPE score, the model accuracy significantly improved from adding the polynomial trend to the regression model, as the error on the testing set dropped from 9.2% to 4.5%:
mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2[1] 0.03688770 0.0421261840.2.4 The tslm Function
So far, we have seen the manual process of transforming a ts object to a linear regression forecasting model format. The tslmfunction from the forecast package provides a built-in function for transforming a ts object into a linear regression forecasting model. Using the tslm function, you can set the regression component along with other features.
We will now repeat the previous example and forecast the last 12 observations of the \(USgas\) series with the tslm function using the trend, square of the trend, and the seasonal components. First, let’s split the series to training and testing partitions using the ts_split function:
Example:
USgas_split <- ts_split(USgas, sample.out = h)
train.ts <- USgas_split$train
test.ts <- USgas_split$testNext, we will apply the same formula we used to create the preceding \(md2\) forecasting model using the tslm function:
library(forecast)
md3 <- tslm(train.ts ~ season + trend + I(trend^2))Let’s now review \(md3\), the output of the tslm function, and compare it with the output of \(md2\):
summary(md3)
Call:
tslm(formula = train.ts ~ season + trend + I(trend^2))
Residuals:
Min 1Q Median 3Q Max
-468.47 -54.66 -2.21 63.11 294.32
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.635e+03 3.224e+01 81.738 < 2e-16 ***
season2 -3.004e+02 3.540e+01 -8.487 3.69e-15 ***
season3 -4.841e+02 3.540e+01 -13.676 < 2e-16 ***
season4 -9.128e+02 3.540e+01 -25.787 < 2e-16 ***
season5 -1.099e+03 3.540e+01 -31.033 < 2e-16 ***
season6 -1.118e+03 3.540e+01 -31.588 < 2e-16 ***
season7 -9.547e+02 3.540e+01 -26.968 < 2e-16 ***
season8 -9.322e+02 3.541e+01 -26.329 < 2e-16 ***
season9 -1.136e+03 3.541e+01 -32.070 < 2e-16 ***
season10 -1.046e+03 3.541e+01 -29.532 < 2e-16 ***
season11 -8.001e+02 3.590e+01 -22.286 < 2e-16 ***
season12 -2.618e+02 3.590e+01 -7.293 5.95e-12 ***
trend -1.270e+00 4.472e-01 -2.840 0.00494 **
I(trend^2) 1.713e-02 1.908e-03 8.977 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 109.1 on 212 degrees of freedom
Multiple R-squared: 0.9379, Adjusted R-squared: 0.9341
F-statistic: 246.1 on 13 and 212 DF, p-value: < 2.2e-16As you can observe from the preceding output, both models (\(md2\) and \(md3\)) are identical.
There are several advantages to using the tslm function, as opposed to manually setting a regression model for the series with the lm function:
Efficiency —does not require transforming the series to a data.frame object and feature engineering.
The output object supports all the functionality of the forecast (such as the
accuracyandcheckresidualsfunctions) and TSstudio packages (such as thetest_forecastandplot_forecastfunctions).
40.2.5 Forecasting a series with multiseasonality components
We will use the UKgrid series, which represents the national grid demand for electricity in the UK and it is available in the UKgrid package, to demonstrate the forecasting approach of a multiseasonality series with a linear regression model.
Example:
library(UKgrid)
# Set the series to daily frequency using the data.frame structure
UKdaily <- extract_grid(type = "data.frame",
columns = "ND",
aggregate = "daily")
class(UKdaily)[1] "data.frame"head(UKdaily) TIMESTAMP ND
1 2005-04-01 1920069
2 2005-04-02 1674699
3 2005-04-03 1631352
4 2005-04-04 1916693
5 2005-04-05 1952082
6 2005-04-06 1964584As you can see, this series has two variables:
\(TIMESTAMP\): A date object used as the series timestamp or index.
\(ND\): The net demand of electricity.
We will use the ts_plot function to plot and review the series structure:
ts_plot(UKdaily,
title = "The UK National Demand for Electricity",
Ytitle = "MW",
Xtitle = "Year")As you can see in the preceding plot, the series has a clear downtrend and has multiple seasonality patterns:
Daily: A cycle of 365 days a year.
Day of the week: A 7-day cycle.
Monthly: Effected from the weather.
Evidence for those patterns can be seen in the following heatmap of the series since 2016 using the ts_heatmap function from the TSstudio package:
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
title = "UK the Daily National Grid Demand Heatmap")As you can see in the series heatmap, the overall demand increases throughout the winter weeks (for example, calendar weeks 1 to 12 and weeks 44 to 52). In addition, you can observe the change of the series during the days of the weeks, as the demand increases during the working days of the week, and decreases over the weekend.
40.3 Alternative Trend Shapes
40.3.1 Exponential Trend
Besides linear trend, several alternative trend shapes can be fitted via a linear regression model. An exponential trend implies a multiplicative increase/decrease of the series over time.
To fit an exponential trend, simply replace the output variable y with \(log(y)\) and fit a linear regression:
\[log(y_t) = \beta_0 + \beta_1t + \epsilon\] Exponential trends are popular in sales data, where they reflect percentage growth.
40.3.2 Polynomial Trend
Another nonlinear trend shape that is easy to fit via linear regression is a polynomial trend, and in particular, a **quadratic relationship of the form:
\[y_t = \beta_0 + \beta_1t + \beta_2t^2 + \epsilon\]
This is done by creating an additional predictor \(t^2\) (the square of \(t\)) and fitting a multiple linear regression with the two predictors \(t\) and \(t^2\).