Time Series Forecasting in R: Methods, Packages, and Datasets
Time series forecasting is an important part of data science, helping us predict future values from observations collected over time. Applications range from forecasting retail sales and electricity demand to analyzing financial markets, monitoring industrial processes, and estimating economic indicators. Unlike traditional machine learning problems, time series forecasting requires us to consider temporal dependencies, trends, seasonal patterns, and changes in the behavior of data.
R provides a comprehensive environment for time series analysis, with specialized packages for statistical modeling, decomposition, visualization, and forecast evaluation. Its built-in datasets and external data collections also make it possible to explore different forecasting techniques without immediately requiring a custom dataset. From classical methods such as exponential smoothing and ARIMA to modern approaches involving regression, decomposition, and machine learning, R supports a wide range of forecasting workflows.
In this guide, we explore Time Series Forecasting in R through its major methods, essential packages, and commonly used datasets. We also examine practical R code examples and explain how to select, train, and evaluate forecasting models for different types of time-dependent data.
Understanding Time Series Forecasting in R
A time series is a sequence of observations recorded at successive points in time. These observations may be collected hourly, daily, weekly, monthly, quarterly, or annually. The objective of forecasting is to use historical observations to estimate future values while accounting for patterns in the underlying data. For example, a retailer might use monthly sales records to estimate demand for the next quarter. An energy provider might forecast hourly electricity consumption to plan resource allocation, while an economist might predict inflation or employment trends using historical economic indicators.
Several characteristics influence the choice of forecasting method:
Trend: A long-term upward or downward movement in the observed values.
Seasonality: A recurring pattern associated with a fixed calendar period, such as increased retail sales during particular months.
Cyclic behavior: Longer-term fluctuations that do not necessarily follow a fixed seasonal schedule.
Autocorrelation: Dependence between observations recorded at different points in time.
Stationarity: A property of a time series whose statistical characteristics, such as its mean and variance, remain stable over time.
Irregular variation: Random fluctuations that remain after the more systematic patterns have been considered.
Understanding these characteristics helps us determine the appropriate forecasting model. For example, a series with strong monthly seasonality may benefit from seasonal ARIMA or exponential smoothing, while a relatively stable series might be adequately forecast using a simple baseline method. R's built-in ts() function provides a convenient way to represent regularly spaced time series data. The function records the starting period and observation frequency, allowing many statistical methods to recognize the temporal structure of the data.
The following example uses the built-in AirPassengers dataset, which contains monthly international airline passenger totals from 1949 through 1960. It is particularly useful for studying trend, seasonality, and changing variation over time.
# Load the built-in airline passenger dataset
data("AirPassengers")
# Display the time series
print(AirPassengers)
# Examine its structure
start(AirPassengers)
end(AirPassengers)
frequency(AirPassengers)
# Visualize the series
plot(
AirPassengers,
main = "Monthly International Airline Passengers",
xlab = "Year",
ylab = "Passengers (thousands)",
col = "steelblue",
lwd = 2
)
# Decompose the series into its components
decomposed <- decompose(AirPassengers)
plot(decomposed)
Output:
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 ...The AirPassengers series has a frequency of 12 because it contains monthly observations. Its decomposition separates the observed data into trend, seasonal, and random components.
The increasing amplitude of its seasonal fluctuations also suggests that a multiplicative decomposition or a logarithmic transformation may be appropriate. Exploratory analysis should generally precede model fitting. Plotting the observations, examining seasonal patterns, and investigating missing values can reveal characteristics that affect model selection. For more complex series, autocorrelation plots and statistical tests can provide additional information about temporal dependence and stationarity.
Major Time Series Forecasting Methods and R Packages
R supports several established forecasting methods, each designed to capture different characteristics of time-dependent data. Rather than selecting the most complicated algorithm, we should begin with methods that provide a reasonable baseline and introduce additional complexity when the data supports it.
1. Baseline forecasting methods
Simple forecasting methods are useful for establishing a reference point against which more sophisticated models can be evaluated. A mean forecast predicts future observations using the historical average. A naïve forecast assumes that the latest observation will remain unchanged, while a seasonal naïve forecast repeats the value from the corresponding season in the previous cycle. The drift method extends the average historical change into the future. These approaches are computationally inexpensive and can perform surprisingly well on suitable datasets. They also help us determine if a complex model delivers a meaningful improvement over a straightforward alternative.
The forecast package provides functions for these methods, alongside more advanced forecasting techniques. The following example creates forecasts from the airline passenger series.
library(forecast)
# Mean forecast
mean_model <- meanf(AirPassengers, h = 12)
# Naive forecast
naive_model <- naive(AirPassengers, h = 12)
# Seasonal naive forecast
seasonal_model <- snaive(AirPassengers, h = 12)
# Drift forecast
drift_model <- rwf(
AirPassengers,
h = 12,
drift = TRUE
)
# Compare the forecast plots
autoplot(seasonal_model)The argument h = 12 requests forecasts for the next 12 monthly observations. For the airline passenger series, seasonal naïve forecasting provides a useful baseline because passenger demand follows a recurring annual pattern. However, its forecasts do not independently estimate the continuing growth in passenger numbers.
Baseline methods are particularly valuable when forecasting systems must be tested regularly. If a complex model cannot consistently outperform an appropriate baseline on unseen observations, its additional complexity may not be justified.
2. Exponential smoothing and ETS models
Exponential smoothing methods generate forecasts using weighted combinations of historical observations, assigning greater importance to recent values in many configurations. ETS models extend this framework by modeling error, trend, and seasonality components. The forecast package supports exponential smoothing through the ets() function. The fable ecosystem provides a modern alternative for modeling tidy time series data. ETS models are useful when the series contains identifiable trends or seasonal patterns that can be represented through level, trend, and seasonal components. They are widely applied to demand forecasting, inventory planning, and operational metrics.
library(forecast)
# Fit an automatically selected ETS model
ets_model <- ets(AirPassengers)
# Generate a 12-month forecast
ets_forecast <- forecast(
ets_model,
h = 12
)
# Display the forecasts and prediction intervals
autoplot(ets_forecast)
# Inspect the fitted model
summary(ets_model)Lets break down the output of this code step by step. The ets() function automatically selects an exponential smoothing model based on the structure of the AirPassengers time series. The fitted model is represented as:
ETS(M,Ad,M)This notation describes three components: multiplicative errors (M), an additive damped trend (Ad), and multiplicative seasonality (M). The model accounts for seasonal fluctuations that change proportionally with the series level while allowing the upward trend to gradually flatten over time. This combination is suitable for the airline passenger dataset, which exhibits both long-term growth and increasing seasonal variation.
The next part of the output from the code above presents the estimated smoothing parameters:
Smoothing parameters:
alpha = 0.7096
beta = 0.0204
gamma = 1e-04
phi = 0.98The smoothing parameters control how the model updates its components as new observations become available. Alpha (0.7096) governs the level and gives recent observations considerable influence. Beta (0.0204) controls trend updates, while gamma (0.0001) indicates very limited updating of the seasonal component. The damping parameter phi (0.98) controls how the trend's influence diminishes over longer forecast horizons.
The model also reports its initial states, which provide the starting values for the level, trend, and seasonal components.
Initial states:
l = 120.9939
b = 1.7705
s = 0.8944 0.7993 0.9217 1.0592 1.2203 1.2318
1.1105 0.9786 0.9804 1.011 0.8869 0.9059
sigma: 0.0392Here, l represents the initial level, and b represents the initial trend. The 12 seasonal values correspond to the monthly seasonal factors. Factors above 1 indicate months with observations above the estimated underlying level, while factors below 1 indicate relatively lower observations. The sigma value represents the estimated scale of the multiplicative error component.
Next, the output provides information criteria that can help compare the fitted model with alternative ETS specifications.
AIC AICc BIC
1395.166 1400.638 1448.623AIC, AICc, and BIC balance model fit against model complexity. Lower values generally indicate a preferable model when comparing candidates fitted to the same dataset using comparable likelihood calculations. AICc applies a correction for sample size, while BIC generally imposes a stronger penalty for model complexity. These values are most useful for model comparison rather than as standalone measures of forecast accuracy.
Finally, the training set error measures summarize how closely the fitted model reproduces the historical observations.
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.5674 10.7473 7.7916 0.4358 2.8579 0.2433 0.0395The mean absolute percentage error (MAPE) is approximately 2.86%, indicating relatively small average absolute percentage errors on the training data. The mean absolute error (MAE) is approximately 7.79, and the root mean squared error (RMSE) is around 10.75, measured in the dataset's original units of thousands of passengers. RMSE penalizes larger errors more heavily than MAE. The mean error (ME) is positive, indicating a small average underforecasting tendency under the usual actual-minus-fitted error convention. The MASE value of 0.2433 indicates that the scaled training errors are below the corresponding naïve benchmark scale. Meanwhile, ACF1 is close to zero, suggesting little residual autocorrelation at lag one.
These training metrics provide a useful initial assessment, but they do not establish how accurately the model will predict future observations. To evaluate generalization, we should test the model on a held-out time period or use rolling-origin cross-validation.
3. ARIMA and seasonal ARIMA
Autoregressive Integrated Moving Average, or ARIMA, is one of the most widely used families of statistical forecasting models. It combines autoregressive terms, differencing, and moving-average terms to represent temporal dependencies in a series.
An ARIMA model is commonly expressed as ARIMA(p, d, q), where:
p represents the autoregressive order.
d represents the degree of differencing.
q represents the moving-average order.
Seasonal ARIMA, often written as SARIMA, extends this structure by incorporating seasonal autoregressive, differencing, and moving-average terms. It is especially useful for data with recurring patterns, such as monthly sales, quarterly economic indicators, or daily demand with weekly seasonality.
The forecast package offers auto.arima() for automatic model selection, while fable provides ARIMA modeling within a tidy workflow. The tseries package also supplies statistical tools, including the Augmented Dickey-Fuller test, that can help investigate stationarity.
library(forecast)
# Automatically select an ARIMA model
arima_model <- auto.arima(
AirPassengers,
seasonal = TRUE
)
# Forecast the next 12 months
arima_forecast <- forecast(
arima_model,
h = 12
)
# Plot the forecast
autoplot(arima_forecast)
# Inspect model coefficients and diagnostics
summary(arima_model)
checkresiduals(arima_model)The auto.arima() function selects an ARIMA model by evaluating candidate specifications for the AirPassengers time series. In this example, the selected model is:
Series: AirPassengers
ARIMA(2,1,1)(0,1,0)[12]The model contains two non-seasonal autoregressive terms, one order of non-seasonal differencing, and one moving-average term. Its seasonal component applies differencing at lag 12, accounting for the annual pattern in monthly observations. The seasonal autoregressive and moving-average orders are both zero. This specification allows the model to capture short-term temporal relationships while addressing trend and yearly seasonality.
The next section reports the estimated model coefficients and their standard errors:
Coefficients:
ar1 ar2 ma1
0.5960 0.2143 -0.9819
s.e. 0.0888 0.0880 0.0292The coefficients ar1 and ar2 represent the first two autoregressive terms, while ma1 represents the first moving-average term. Their estimates describe how the differenced series and past forecast errors contribute to the model. The standard errors indicate the estimation uncertainty associated with these coefficients. The MA coefficient of -0.9819 is close to the boundary of the usual invertibility region, so it is worth checking the fitted model's diagnostics rather than interpreting the coefficient in isolation.
The output then provides the estimated residual variance and the model-selection statistics:
sigma^2 = 132.3
log likelihood = -504.92
AIC = 1017.85
AICc = 1018.17
BIC = 1029.35The residual variance (sigma² = 132.3) measures the remaining variation captured by the model's error term. The log likelihood describes how well the model fits the observed data under its statistical assumptions. AIC, AICc, and BIC balance model fit against model complexity, with lower values generally preferred when comparing suitable candidate models fitted to the same series. These statistics are useful for model selection but do not directly measure forecast accuracy on unseen observations.
Finally, the training set error measures summarize how closely the fitted model reproduces the historical passenger data:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.3423 10.8462 7.8675 0.4207 2.8005 0.2456 -0.0012The mean absolute percentage error (MAPE) is approximately 2.80%, indicating a relatively small average absolute percentage error on the training data. The mean absolute error (MAE) is around 7.87, while the root mean squared error (RMSE) is approximately 10.85, both measured in thousands of passengers. The positive mean error (ME) indicates a small average underforecasting tendency under the usual actual-minus-fitted error convention. The MASE value of 0.2456 indicates that the scaled training error is below its naïve benchmark scale. Finally, ACF1 is close to zero, suggesting little residual autocorrelation at lag one, although additional residual checks are needed to assess the model fully.
Evaluating model performance is a crucial step in time series forecasting, and analyzing residual diagnostics tells us exactly how well a model has captured the underlying data patterns. This image showcases the standard diagnostic plots for a SARIMA(2,1,1)(0,1,0)[12] model evaluated from 1949 to 1961.
The top time-series plot shows the residuals generally hovering around a zero mean, though a noticeable jump in volatility toward 1960 hints at potential heteroscedasticity. Meanwhile, the bottom-left ACF plot checks for any lingering patterns; while most lags stay safely within the blue white-noise boundaries, a prominent spike breaks through near lag 23, signaling that some seasonal correlation has been left on the table.
Looking at the final chart, we can assess our error distribution. The bottom-right histogram maps the residuals against a theoretical normal curve, showing a mostly symmetric, bell-shaped distribution around zero that validates our normality assumption despite a few extreme outliers. Ultimately, these diagnostic visualizations serve as essential check-ins to prove whether our model is truly ready for prime-time forecasting or if it requires a few parameter tweaks to clean up the remaining seasonal patterns.
The Ljung–Box test evaluates the residuals to determine if significant autocorrelation remains after fitting the ARIMA model. The output is shown below.
Ljung-Box test
data: Residuals from ARIMA(2,1,1)(0,1,0)[12]
Q* = 37.784, df = 21, p-value = 0.01366
Model df: 3. Total lags used: 24The test produces a Q-statistic of 37.784 with 21 degrees of freedom and a p-value of 0.01366. The null hypothesis states that the residuals have no significant autocorrelation up to the tested lags. Since the p-value is below the commonly used significance level of 0.05, we reject the null hypothesis. This suggests that the ARIMA model has not captured all the temporal dependencies in the data, even though its lag-one residual autocorrelation is close to zero.
The test uses 24 lags, with the degrees of freedom adjusted by the three estimated AR and MA coefficients. These results indicate that we should investigate the residual autocorrelation further and consider alternative model specifications. Although the fitted model shows reasonable training accuracy, the Ljung–Box test suggests that its residuals are not entirely random, highlighting the importance of checking model diagnostics alongside forecast error measures.
4. Decomposition, regression, and specialized models
Not every forecasting problem is best handled by a single statistical model. Some datasets contain multiple seasonal patterns, external predictors, changing trends, or dependencies across several related variables. R offers several additional approaches:
STL decomposition: Separates a series into seasonal, trend, and remainder components. The stats package provides stl(), while the forecast package includes stlf() for forecasting using STL decomposition.
Dynamic regression: Combines regression with time series error structures. It is useful when external predictors, such as prices, promotions, or weather measurements, influence the target series.
TBATS: Supports complex seasonal structures, including multiple seasonal periods. It is available through the forecast package.
Prophet: Provides trend, seasonality, and holiday components through the prophet package. It can be useful for business series with calendar effects and changing trends.
VAR models: Vector autoregression models multiple related time series and their lagged relationships. The vars package supports this approach.
GARCH models: The rugarch package provides tools for modeling changing conditional variance, particularly in financial time series. These models are often used to forecast volatility rather than the mean alone.
The appropriate choice depends on the data and the forecasting objective. For example, STL can help separate a changing trend from recurring seasonal behavior, while dynamic regression allows us to incorporate explanatory variables that may improve forecasts. TBATS can be useful when multiple seasonal patterns make simpler seasonal models difficult to apply.
The most important principle is to match the model to the structure of the series. A specialized algorithm is not automatically better than a simpler method, and every candidate should be evaluated using observations that were not used during model fitting.
Essential R Packages and Datasets for Forecasting
The strength of R's time series ecosystem comes from the range of packages available for different stages of the forecasting workflow. Some packages focus on classical statistical models, while others provide tidy data structures, visualization tools, specialized algorithms, or collections of datasets for experimentation. The forecast package remains a useful starting point for traditional forecasting methods. It includes functions for ARIMA, ETS, baseline forecasts, decomposition-based forecasting, and forecast accuracy evaluation. Its ts-oriented workflow is convenient for learning classical techniques.
The fpp3 package provides a different approach through the tidyverts ecosystem. It loads several related packages, including tsibble, fable, and feasts, which support data management, forecasting models, and time series analysis. Here are some important packages to consider:
forecast: Classical forecasting methods, automatic ARIMA selection, ETS models, TBATS, and accuracy metrics.
fable: Tidy forecasting models, including ARIMA and ETS, with support for model comparison and multiple time series.
tsibble: Represents time series data in a tidy structure with explicit time indexes.
feasts: Provides time series visualization, decomposition, and feature extraction tools.
tseries: Includes statistical tests such as the Augmented Dickey-Fuller test and other time series utilities.
prophet: Supports trend, seasonality, and holiday-based forecasting.
vars: Implements vector autoregression and related multivariate time series methods.
rugarch: Supports generalized autoregressive conditional heteroskedasticity models for volatility analysis.
For a traditional forecasting workflow, forecast is a practical starting point. For projects involving multiple series, grouped analysis, or tidy data transformations, fable and tsibble offer a structured alternative. The two approaches can complement each other, although their data representations and modeling interfaces differ.
Working with datasets available in R
Built-in datasets make it easier to learn forecasting techniques because they provide ready-to-use examples with established temporal structures. Several are available through R's datasets package, which is included with standard R installations. Others become available after installing specialized forecasting packages. Commonly used datasets include:
AirPassengers: Monthly international airline passenger totals, useful for studying trend and annual seasonality.
nottem: Average air temperatures at Nottingham Castle, recorded monthly from 1920 through 1939.
Nile: Annual measurements of the flow of the Nile River near Aswan, useful for studying long-term changes.
JohnsonJohnson: Quarterly earnings per share for Johnson & Johnson, useful for examining trend and seasonality.
LakeHuron: Annual measurements of the level of Lake Huron, suitable for exploring long-term variation.
aus_retail: Australian retail turnover data available through the fpp3 ecosystem, useful for exploring retail demand and seasonal patterns.
tourism: Australian tourism data available through the same ecosystem, suitable for grouped forecasting and seasonal analysis.
vic_elec: Victorian electricity demand data available through the fpp3 ecosystem, useful for high-frequency forecasting and multiple seasonal patterns.
The built-in datasets can be accessed directly after loading the relevant package. For datasets provided by fpp3, install and load the package first.
The following example demonstrates how to work with the aus_retail dataset using the tidyverts ecosystem.
# Install once if the package is not available
install.packages("fpp3")
library(fpp3)
# Load the retail turnover dataset
data("aus_retail")
# Inspect the dataset
glimpse(aus_retail)
# Select one industry and one state
retail_subset <- aus_retail |>
filter(
State == "New South Wales",
Industry == "Food retailing"
)
# Visualize retail turnover over time
retail_subset |>
autoplot(Turnover) +
labs(
title = "Food Retail Turnover in New South Wales",
x = "Year",
y = "Turnover"
)The aus_retail dataset contains multiple time series identified by variables such as state and industry. Filtering it to one state and industry produces a more focused series for exploratory analysis. Its tidy structure also makes it easier to compare categories or extend the workflow to grouped forecasting.
Before fitting a model, we should verify the observation frequency, identify missing periods, inspect the measurement units, and check for revisions or unusual values. A dataset containing several locations or product categories may require separate models or a global forecasting approach that learns across multiple series.
Dataset selection should reflect the forecasting problem. AirPassengers is useful for introductory seasonal forecasting, Nile is appropriate for annual observations, and electricity demand data provides a more challenging setting involving high-frequency observations and multiple seasonal cycles.
Building, Evaluating, and Improving Forecasting Models in R
A reliable forecasting workflow involves more than fitting a model to historical observations. We must also establish a suitable forecast horizon, separate training and testing periods, evaluate prediction errors, and investigate the model's behavior under realistic forecasting conditions. Randomly splitting time series observations into training and testing sets is generally inappropriate because it can allow future observations to influence the evaluation of past predictions. Instead, we should preserve chronological order. The earlier observations are used for model training, while the later observations are reserved for testing.
For example, when forecasting monthly sales, we might train a model on the first several years and evaluate it on the final 12 months. This setup provides a more realistic estimate of performance when forecasting future months from currently available information.
The forecast package offers several accuracy measures, including mean absolute error (MAE), root mean squared error (RMSE), and mean absolute percentage error (MAPE). Each captures a different aspect of forecasting performance.
MAE: Measures the average absolute difference between actual and forecast values. It is expressed in the same units as the target variable.
RMSE: Penalizes larger errors more heavily because the errors are squared before averaging.
MAPE: Expresses absolute errors as percentages of actual values, but it becomes undefined or unstable when actual observations are zero or close to zero.
MASE: Scales forecast errors against a naïve benchmark and can support comparisons across series with different units.
No single metric is suitable for every forecasting task. When large errors are especially costly, RMSE may be informative. When straightforward interpretation in the target's original units is important, MAE is often useful. For intermittent demand or data containing zeros, MAPE should be avoided or supplemented with other measures.
The following example demonstrates a chronological train-test split using the airline passenger dataset. It fits an ETS model on the training period and evaluates forecasts against the final 24 months.
library(forecast)
# Reserve the final 24 months for testing
train <- head(
AirPassengers,
length(AirPassengers) - 24
)
test <- tail(
AirPassengers,
24
)
# Fit the model using training observations only
model <- ets(train)
# Forecast the test period
predictions <- forecast(
model,
h = length(test)
)
# Compare forecasts with actual observations
accuracy(predictions, test)
# Plot forecasts and prediction intervals
autoplot(predictions) +
autolayer(
test,
series = "Actual",
color = "red"
)
Output:
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set 1.208427 8.898232 6.65393 0.3935913 2.777234 0.232866 0.1188263 NA
Test set 63.210920 72.547909 63.21297 13.3028462 13.303447 2.212249 0.7463877 1.356753This evaluation estimates how well the fitted model forecasts the final two years of the dataset. However, a single test period may not represent performance across different market conditions or seasonal cycles. For more robust comparisons, rolling-origin evaluation repeatedly moves the forecast origin forward, refits the model using the observations available at each point, and evaluates subsequent forecasts.
After comparing several candidate models, we should also inspect residuals and consider the operational cost of maintaining each approach. A model that performs slightly better but requires extensive data preparation or frequent manual intervention may not be the best practical choice.
The upper red boundary marks the outer edge of the 95% prediction interval, highlighting the model's maximum projected variance. Because the model uses multiplicative seasonality, the forecasted seasonal peaks expand dramatically relative to the upward trend, creating a wide gap between the average forecast (blue line) and the upper threshold. This red line represents the statistical worst-case scenario for high demand, visually warning stakeholders of potential extreme spikes in the forecast horizon.
For production forecasting, the workflow should include data validation, scheduled model fitting, monitoring of forecast errors, and periodic reassessment of model assumptions. Forecasts should also be compared with simple baselines to detect performance deterioration.
Finally, forecasting results should be communicated with appropriate uncertainty estimates. Prediction intervals help decision-makers understand the range of plausible future outcomes rather than treating a point forecast as a guaranteed result. Their reliability depends on the fitted model, the data-generating process, and the assumption that historical patterns remain informative about the future.
Conclusion
Time Series Forecasting in R combines statistical methods, specialized packages, and accessible datasets into a flexible framework for predicting future observations. Baseline methods establish performance benchmarks, ETS and ARIMA models capture important temporal structures, and approaches such as STL, dynamic regression, TBATS, Prophet, and VAR address more specialized forecasting requirements.
R's built-in datasets provide a convenient starting point for experimentation, while the forecast and tidyverts ecosystems support increasingly complex analytical workflows. By understanding the characteristics of the data, selecting appropriate methods, and evaluating forecasts on unseen observations, we can develop forecasting systems that are both technically sound and useful in real-world applications.
If you want to fully master time series analysis and forecasting, our comprehensive reference guide explores the entire landscape of temporal modeling. Moving beyond just code, it breaks down essential forecasting methodologies including univariate, multivariate, and intermittent structures alongside standard industry benchmark datasets spanning retail, finance, and energy systems. To dive deep into the foundational concepts, models, and practical frameworks of chronological data, check out the full article here: Time Series Forecasting: Models, Techniques, and Examples.
