ARCH Model
Stepbystep guide on using the ARCH Model
with Statsforecast
.
Table of Contents
 Introduction
 ARCH Models
 Loading libraries and data
 Explore data with the plot method
 Split the data into training and testing
 Implementation of ARCH with StatsForecast
 Crossvalidation
 Model evaluation
 References
Introduction
Financial time series analysis has been one of the hottest research
topics in the recent decades. In this guide, we illustrate the stylized
facts of financial time series by real financial data. To characterize
these facts, new models different from the Box Jenkins ones are needed.
And for this reason, ARCH models were firstly proposed by R. F. Engle in
1982 and have been extended by a great number of scholars since then. We
also demonstrate how to use Python and its libraries to implement
ARCH
.
As we have known, there are lot of time series that possess the ARCH effect, that is, although the (modeling residual) series is white noise, its squared series may be autocorrelated. What is more, in practice, a large number of financial time series are found having this property so that the ARCH effect has become one of the stylized facts from financial time series.
Stylized Facts of Financial Time Series
Now we briefly list and describe several important stylized facts (features) of financial return series:

Fat (heavy) tails: The distribution density function of returns often has fatter (heavier) tails than the tails of the corresponding normal distribution density.

ARCH effect: Although the return series can often be seen as a white noise, its squared (and absolute) series may usually be autocorrelated, and these autocorrelations are hardly negative.

Volatility clustering: Large changes in returns tend to cluster in time, and small changes tend to be followed by small changes.

Asymmetry: As we have know , the distribution of asset returns is slightly negatively skewed. One possible explanation could be that traders react more strongly to unfavorable information than favorable information.
Definition of ARCH Models
Specifically, we give the definition of the ARCH model as follows.
Definition 1. An $\text{ARCH(p)}$ model with order $p≥1$ is of the form
where $\omega ≥ 0, \alpha_i ≥ 0$, and $\alpha_p > 0$ are constants, $\varepsilon_t \sim iid(0, 1)$, and $\varepsilon_t$ is independent of $\{X_k;k ≤ t − 1 \}$. A stochastic process $X_t$ is called an $ARCH(p)$ process if it satisfies Eq. (1).
By Definition 1, $\sigma_{t}^2$ (and $\sigma_t$ ) is independent of $\varepsilon_t$ . Besides, usually it is further assumed that $\varepsilon_t \sim N(0, 1)$. Sometimes, however, we need to further suppose that $\varepsilon_t$ follows a standardized (skew) Student’s T distribution or a generalized error distribution in order to capture more features of a financial time series.
Let $\mathscr{F}_s$ denote the information set generated by $\{X_k;k ≤ s \}$, namely, the sigma field $\sigma(X_k;k ≤ s)$. It is easy to see that $\mathscr{F}_s$ is independent of $\varepsilon_t$ for any $s <t$. According to Definition 1 and the properties of the conditional mathematical expectation, we have that
and
This implies that $\sigma_{t}^2$ is the conditional variance of $X_t$ and it evolves according to the previous values of $\{X_{k}^2; t −p ≤ k ≤ t −1\}$ like an $\text{AR}(p)$ model. And so Model (1) is named an $\text{ARCH}(p)$ model.
As an example of $\text{ARCH}(p)$ models, let us consider the $\text{ARCH(1)}$ model
Explicitly, the unconditional mean $E(X_t) = E(\sigma_t \varepsilon_t) = E(\sigma_t) E(\varepsilon_t) = 0.$
Additionally, the ARCH(1) model can be expressed as
$X_{t}^2 =\sigma_{t}^2 +X_{t}^2 − \sigma_{t}^2 =\omega +\alpha_1 X_{t1}^2 +\sigma_{t}^2 \varepsilon_{t}^2 −\sigma_{t}^2 =\omega +\alpha_1 X_{t}^2 +\eta_t$
that is,
where $\eta_t = \sigma_{t}^2(\varepsilon_{t}^2 − 1)$. It can been shown that $\eta_t$ is a new white noise, which is left as an exercise for reader. Hence, if $0 < \alpha_1 < 1$, Eq. (4) is a stationary $\text{AR(1)}$ model for the series Xt2. Thus, the unconditional variance
$Var ( X_t ) = E( X_{t}^2 ) = E(\omega+ \alpha_1 X_{t1}^2 + \eta_t ) = \omega+ \alpha_1 E( X_{t}^2 ) ,$
that is, $Var (X_t) = E (X_{t}^2 ) =\frac{\omega}{1\alpha_1}$
Moreover, for $h > 0$, in light of the properties of the conditional mathematical expectation and by (2), we have that
$E(X_{t+h} X_t) = E(E(X_{t+h} X_t\mathscr{F}_{t+h1})) = E(X_t E(X_{t+h}\mathscr{F}_{t+h1})) = 0.$
In conclusion, if $0 < \alpha_1 < 1$, we have that:

Any $\text{ARCH}(1)$ process $\{X_t \}$ defined by Eqs.(3) follows a white noise $WN(0, \omega/(1 − \alpha_1))$ .

Since $X_{t}^2$ is an $\text{AR}(1)$ process defined by (4), $\text{Corr}(X_{t}^2,X_{t+h}^2) = \alpha_{1}^{h} > 0$, which reveals the ARCH effect.

It is clear that $E(\eta_t\mathscr{F}_s)=0$ for any $t>s$,and with Eq.(4),for any $k>1$: $Var(X_{t+k} \mathscr{F}_t ) = E(X_{t+K}^2 \mathscr{F}_t)$ $= E(\omega + \alpha_1 X_{t+k1}+ \eta_{t+k}\mathscr{F}_t )$ $= \omega + \alpha_1 Var(X_{t+k−1}\mathscr{F}_t),$
which reflects the volatility clustering, that is, large (small) volatility is followed by large (small) one.
In addition, we are able to prove that Xt defined by Eq. (3) has heavier tails than the corresponding normal distribution. At last, note that these properties of the ARCH(1) model can be generalized to ARCH(p) models.
Advantages and disadvantages of the Autoregressive Conditional Heteroskedasticity (ARCH) model:
Advantages  Disadvantages 

 The ARCH model is useful for modeling volatility in financial time series, which is important for investment decision making and risk management.   The ARCH model assumes that the forecast errors are independent and identically distributed, which may not be realistic in some cases. 
 The ARCH model takes heteroscedasticity into account, which means that it can model time series with variances that change over time.   The ARCH model can be difficult to fit to data with many parameters, which may require large amounts of data or advanced estimation techniques. 
 The ARCH model is relatively easy to use and can be implemented with standard econometrics software.   The ARCH model does not take into account the possible relationship between the mean and the variance of the time series, which may be important in some cases. 
Note:
The ARCH model is a useful tool for modeling volatility in financial time series, but like any econometric model, it has limitations and should be used with caution depending on the specific characteristics of the data being modeled.
Autoregressive Conditional Heteroskedasticity (ARCH) Applications

Finance  The ARCH model is widely used in finance to model volatility in financial time series, such as stock prices, exchange rates, interest rates, etc.

Economics  The ARCH model can be used to model volatility in economic data, such as GDP, inflation, unemployment, among others.

Engineering  The ARCH model can be used in engineering to model volatility in data related to energy, climate, pollution, industrial production, among others.

Social Sciences  The ARCH model can be used in the social sciences to model volatility in data related to demography, health, education, among others.

Biology  The ARCH model can be used in biology to model volatility in data related to evolution, genetics, epidemiology, among others.
Loading libraries and data
Tip
Statsforecast will be needed. To install, see instructions.
Next, we import plotting libraries and configure the plotting style.
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
'figure.facecolor': '#212946',
'axes.facecolor': '#212946',
'savefig.facecolor':'#212946',
'axes.grid': True,
'axes.grid.which': 'both',
'axes.spines.left': False,
'axes.spines.right': False,
'axes.spines.top': False,
'axes.spines.bottom': False,
'grid.color': '#2A3459',
'grid.linewidth': '1',
'text.color': '0.9',
'axes.labelcolor': '0.9',
'xtick.color': '0.9',
'ytick.color': '0.9',
'font.size': 12 }
plt.rcParams.update(dark_style)
from pylab import rcParams
rcParams['figure.figsize'] = (18,7)
Read Data
Let’s pull the S&P500 stock data from the Yahoo Finance site.
import pandas as pd
import time
from datetime import datetime
import yfinance as yf
ticker = '^GSPC'
period1 = int(time.mktime(datetime(2015, 1, 1, 23, 59).timetuple()))
period2 = int(time.mktime(datetime.now().timetuple()))
interval = '1d' # 1d, 1m
SP_500 = yf.download(ticker, start=period1, end=period2, interval=interval)
SP_500 = SP_500.reset_index()
SP_500.head()
[*********************100%***********************] 1 of 1 completed
Date  Open  High  Low  Close  Adj Close  Volume  

0  20150102  2058.899902  2072.360107  2046.040039  2058.199951  2058.199951  2708700000 
1  20150105  2054.439941  2054.439941  2017.339966  2020.579956  2020.579956  3799120000 
2  20150106  2022.150024  2030.250000  1992.439941  2002.609985  2002.609985  4460110000 
3  20150107  2005.550049  2029.609985  2005.550049  2025.900024  2025.900024  3805480000 
4  20150108  2030.609985  2064.080078  2030.609985  2062.139893  2062.139893  3934010000 
df=SP_500[["Date","Close"]]
The input to StatsForecast is always a data frame in long format with three columns: unique_id, ds and y:

The
unique_id
(string, int or category) represents an identifier for the series. 
The
ds
(datestamp) column should be of a format expected by Pandas, ideally YYYYMMDD for a date or YYYYMMDD HH:MM:SS for a timestamp. 
The
y
(numeric) represents the measurement we wish to forecast.
df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()
ds  y  unique_id  

0  20150102  2058.199951  1 
1  20150105  2020.579956  1 
2  20150106  2002.609985  1 
3  20150107  2025.900024  1 
4  20150108  2062.139893  1 
print(df.dtypes)
ds datetime64[ns]
y float64
unique_id object
dtype: object
Explore data with the plot method
Plot a series using the plot method from the StatsForecast class. This method prints a random series from the dataset and is useful for basic EDA.
from statsforecast import StatsForecast
StatsForecast.plot(df)
The Augmented DickeyFuller Test
An Augmented DickeyFuller (ADF) test is a type of statistical test that determines whether a unit root is present in time series data. Unit roots can cause unpredictable results in time series analysis. A null hypothesis is formed in the unit root test to determine how strongly time series data is affected by a trend. By accepting the null hypothesis, we accept the evidence that the time series data is not stationary. By rejecting the null hypothesis or accepting the alternative hypothesis, we accept the evidence that the time series data is generated by a stationary process. This process is also known as stationary trend. The values of the ADF test statistic are negative. Lower ADF values indicate a stronger rejection of the null hypothesis.
Augmented DickeyFuller Test is a common statistical test used to test whether a given time series is stationary or not. We can achieve this by defining the null and alternate hypothesis.
Null Hypothesis: Time Series is nonstationary. It gives a timedependent trend. Alternate Hypothesis: Time Series is stationary. In another term, the series doesn’t depend on time.
ADF or t Statistic < critical values: Reject the null hypothesis, time series is stationary. ADF or t Statistic > critical values: Failed to reject the null hypothesis, time series is nonstationary.
Let’s check if our series that we are analyzing is a stationary series.
Let’s create a function to check, using the Dickey Fuller
test
from statsmodels.tsa.stattools import adfuller
def Augmented_Dickey_Fuller_Test_func(series , column_name):
print (f'DickeyFuller test results for columns: {column_name}')
dftest = adfuller(series, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','pvalue','No Lags Used','Number of observations used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print (dfoutput)
if dftest[1] <= 0.05:
print("Conclusion:====>")
print("Reject the null hypothesis")
print("The data is stationary")
else:
print("Conclusion:====>")
print("The null hypothesis cannot be rejected")
print("The data is not stationary")
Augmented_Dickey_Fuller_Test_func(df["y"],'S&P500')
DickeyFuller test results for columns: S&P500
Test Statistic 0.814971
pvalue 0.814685
No Lags Used 10.000000
...
Critical Value (1%) 3.433341
Critical Value (5%) 2.862861
Critical Value (10%) 2.567473
Length: 7, dtype: float64
Conclusion:====>
The null hypothesis cannot be rejected
The data is not stationary
In the previous result we can see that the Augmented_Dickey_Fuller
test gives us a pvalue
of 0.864700, which tells us that the null
hypothesis cannot be rejected, and on the other hand the data of our
series are not stationary.
We need to differentiate our time series, in order to convert the data to stationary.
Return Series
Since the 1970s, the financial industry has been very prosperous with advancement of computer and Internet technology. Trade of financial products (including various derivatives) generates a huge amount of data which form financial time series. For finance, the return on a financial product is most interesting, and so our attention focuses on the return series. If $P_t$ is the closing price at time t for a certain financial product, then the return on this product is
$X_t = \frac{(P_t − P_{t−1})}{P_{t−1}} ≈ log(P_t ) − log(P_{t−1}).$
It is return series $\{X_t \}$ that have been much independently studied. And important stylized features which are common across many instruments, markets, and time periods have been summarized. Note that if you purchase the financial product, then it becomes your asset, and its returns become your asset returns. Now let us look at the following examples.
We can estimate the series of returns using the
pandas,
DataFrame.pct_change()
function. The pct_change()
function has a
periods parameter whose default value is 1. If you want to calculate a
30day return, you must change the value to 30.
df['return'] = 100 * df["y"].pct_change()
df.dropna(inplace=True, how='any')
df.head()
ds  y  unique_id  return  

1  20150105  2020.579956  1  1.827811 
2  20150106  2002.609985  1  0.889347 
3  20150107  2025.900024  1  1.162984 
4  20150108  2062.139893  1  1.788828 
5  20150109  2044.810059  1  0.840381 
import plotly.express as px
fig = px.line(df, x=df["ds"], y="return",title="SP500 Return Chart",template = "plotly_dark")
fig.show()
Creating Squared Returns
df['sq_return'] = df["return"].mul(df["return"])
df.head()
ds  y  unique_id  return  sq_return  

1  20150105  2020.579956  1  1.827811  3.340891 
2  20150106  2002.609985  1  0.889347  0.790938 
3  20150107  2025.900024  1  1.162984  1.352532 
4  20150108  2062.139893  1  1.788828  3.199906 
5  20150109  2044.810059  1  0.840381  0.706240 
Returns vs Squared Returns
from plotly.subplots import make_subplots
import plotly.graph_objects as go
fig = make_subplots(rows=1, cols=2)
fig.add_trace(go.Scatter(x=df["ds"], y=df["return"],
mode='lines',
name='return'),
row=1, col=1
)
fig.add_trace(go.Scatter(x=df["ds"], y=df["sq_return"],
mode='lines',
name='sq_return'),
row=1, col=2
)
fig.update_layout(height=600, width=800, title_text="Returns vs Squared Returns", template = "plotly_dark")
fig.show()
from scipy.stats import probplot, moment
from statsmodels.tsa.stattools import adfuller, q_stat, acf
import numpy as np
import seaborn as sns
def plot_correlogram(x, lags=None, title=None):
lags = min(10, int(len(x)/5)) if lags is None else lags
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 8))
x.plot(ax=axes[0][0], title='Return')
x.rolling(21).mean().plot(ax=axes[0][0], c='k', lw=1)
q_p = np.max(q_stat(acf(x, nlags=lags), len(x))[1])
stats = f'QStat: {np.max(q_p):>8.2f}\nADF: {adfuller(x)[1]:>11.2f}'
axes[0][0].text(x=.02, y=.85, s=stats, transform=axes[0][0].transAxes)
probplot(x, plot=axes[0][1])
mean, var, skew, kurtosis = moment(x, moment=[1, 2, 3, 4])
s = f'Mean: {mean:>12.2f}\nSD: {np.sqrt(var):>16.2f}\nSkew: {skew:12.2f}\nKurtosis:{kurtosis:9.2f}'
axes[0][1].text(x=.02, y=.75, s=s, transform=axes[0][1].transAxes)
plot_acf(x=x, lags=lags, zero=False, ax=axes[1][0])
plot_pacf(x, lags=lags, zero=False, ax=axes[1][1])
axes[1][0].set_xlabel('Lag')
axes[1][1].set_xlabel('Lag')
fig.suptitle(title+ f'DickeyFuller: {adfuller(x)[1]:>11.2f}', fontsize=14)
sns.despine()
fig.tight_layout()
fig.subplots_adjust(top=.9)
plot_correlogram(df["return"], lags=30, title="Time Series Analysis plot \n")
LjungBox Test
LjungBox is a test for autocorrelation that we can use in tandem with our ACF and PACF plots. The LjungBox test takes our data, optionally either lag values to test, or the largest lag value to consider, and whether to compute the BoxPierce statistic. LjungBox and BoxPierce are two similar test statisitcs, Q , that are compared against a chisquared distribution to determine if the series is white noise. We might use the LjungBox test on the residuals of our model to look for autocorrelation, ideally our residuals would be white noise.
 Ho : The data are independently distributed, no autocorrelation.
 Ha : The data are not independently distributed; they exhibit serial correlation.
The LjungBox with the BoxPierce option will return, for each lag, the LjungBox test statistic, LjungBox pvalues, BoxPierce test statistic, and BoxPierce pvalues.
If $p<\alpha (0.05)$ we reject the null hypothesis.
from statsmodels.stats.diagnostic import acorr_ljungbox
ljung_res = acorr_ljungbox(df["return"], lags= 40, boxpierce=True)
ljung_res.head()
lb_stat  lb_pvalue  bp_stat  bp_pvalue  

1  49.222273  2.285409e12  49.155183  2.364927e12 
2  62.991348  2.097020e14  62.899234  2.195861e14 
3  63.944944  8.433622e14  63.850663  8.834380e14 
4  74.343652  2.742989e15  74.221024  2.911751e15 
5  80.234862  7.494100e16  80.093498  8.022242e16 
Split the data into training and testing
Let’s divide our data into sets
 Data to train our
ARCH
model  Data to test our model
For the test data we will use the last 30 day to test and evaluate the performance of our model.
df=df[["ds","unique_id","return"]]
df.columns=["ds", "unique_id", "y"]
train = df[df.ds<='20230524'] # Let's forecast the last 30 days
test = df[df.ds>'20230524']
train.shape, test.shape
((2112, 3), (87, 3))
Now let’s plot the training data and the test data.
sns.lineplot(train,x="ds", y="y", label="Train")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.show()
Implementation of ARCH with StatsForecast
To also know more about the parameters of the functions of the
ARCH Model
, they are listed below. For more information, visit the
documentation
p : int
Number of lagged versions of the series.
alias : str
Custom name of the model.
prediction_intervals : Optional[ConformalIntervals]
Information to compute conformal prediction intervals.
By default, the model will compute the native prediction
intervals.
Load libraries
from statsforecast import StatsForecast
from statsforecast.models import ARCH
Building Model
Import and instantiate the models. Setting the argument is sometimes tricky. This article on Seasonal periods by the master, Rob Hyndmann, can be useful.season_length.
season_length = 7 # Daily data
horizon = len(test) # number of predictions biasadj=True, include_drift=True,
models = [ARCH(p=2,)]
We fit the models by instantiating a new StatsForecast object with the following parameters:
models: a list of models. Select the models you want from models and import them.

freq:
a string indicating the frequency of the data. (See pandas’ available frequencies.) 
n_jobs:
n_jobs: int, number of jobs used in the parallel processing, use 1 for all cores. 
fallback_model:
a model to be used if a model fails.
Any settings are passed into the constructor. Then you call its fit method and pass in the historical data frame.
sf = StatsForecast(df=train,
models=models,
freq='C', # custom business day frequency
n_jobs=1)
Fit the Model
sf.fit()
StatsForecast(models=[ARCH(2)])
Let’s see the results of our ARCH model. We can observe it with the following instruction:
result=sf.fitted_[0,0].model_
result
{'p': 2,
'q': 0,
'coeff': array([0.44320919, 0.34706759, 0.35171967]),
'message': 'Optimization terminated successfully',
'y_vals': array([1.12220268, 0.73186004]),
'sigma2_vals': array([1.38768694, nan, 1.89277546, ..., 0.76423015, 0.45064543,
0.88036943]),
'fitted': array([ nan, nan, 2.23474473, ..., 1.48032981,
1.10018826, 0.98050094]),
'actual_residuals': array([ nan, nan, 1.07176046, ..., 1.4958333 ,
2.22239094, 0.2486409 ])}
Let us now visualize the residuals of our models.
As we can see, the result obtained above has an output in a dictionary,
to extract each element from the dictionary we are going to use the
.get()
function to extract the element and then we are going to save
it in a pd.DataFrame()
.
residual=pd.DataFrame(result.get("actual_residuals"), columns=["residual Model"])
residual
residual Model  

0  NaN 
1  NaN 
2  1.071760 
…  … 
2109  1.495833 
2110  2.222391 
2111  0.248641 
import scipy.stats as stats
fig, axs = plt.subplots(nrows=2, ncols=2)
# plot[1,1]
residual.plot(ax=axs[0,0])
axs[0,0].set_title("Residuals");
# plot
sns.distplot(residual, ax=axs[0,1]);
axs[0,1].set_title("Density plot  Residual");
# plot
stats.probplot(residual["residual Model"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot QQ')
# plot
plot_acf(residual, lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");
plt.show();
Forecast Method
If you want to gain speed in productive settings where you have multiple
series or models we recommend using the
StatsForecast.forecast
method instead of .fit
and .predict
.
The main difference is that the .forecast
doest not store the fitted
values and is highly scalable in distributed environments.
The forecast method takes two arguments: forecasts next h
(horizon)
and level
.

h (int):
represents the forecast h steps into the future. In this case, 12 months ahead. 
level (list of floats):
this optional parameter is used for probabilistic forecasting. Set the level (or confidence percentile) of your prediction interval. For example,level=[90]
means that the model expects the real value to be inside that interval 90% of the times.
The forecast object here is a new data frame that includes a column with the name of the model and the y hat values, as well as columns for the uncertainty intervals. Depending on your computer, this step should take around 1min.
Y_hat = sf.forecast(horizon, fitted=True)
Y_hat
ds  ARCH(2)  

unique_id  
1  20230525  1.681836 
1  20230526  0.777028 
1  20230529  0.677960 
…  …  … 
1  20230920  0.136752 
1  20230921  0.082173 
1  20230922  0.450958 
values=sf.forecast_fitted_values()
values.head()
ds  y  ARCH(2)  

unique_id  
1  20150105  1.827811  NaN 
1  20150106  0.889347  NaN 
1  20150107  1.162984  2.234745 
1  20150108  1.788828  0.667577 
1  20150109  0.840381  0.752437 
Adding 95% confidence interval with the forecast method
sf.forecast(h=horizon, level=[95])
ds  ARCH(2)  ARCH(2)lo95  ARCH(2)hi95  

unique_id  
1  20230525  1.681836  0.419322  3.782995 
1  20230526  0.777028  3.939044  2.384989 
1  20230529  0.677960  3.907244  2.551323 
…  …  …  …  … 
1  20230920  0.136752  0.795371  1.068876 
1  20230921  0.082173  0.852268  1.016615 
1  20230922  0.450958  1.337117  0.435202 
Y_hat=Y_hat.reset_index()
Y_hat
unique_id  ds  ARCH(2)  

0  1  20230525  1.681836 
1  1  20230526  0.777028 
2  1  20230529  0.677960 
…  …  …  … 
84  1  20230920  0.136752 
85  1  20230921  0.082173 
86  1  20230922  0.450958 
# Merge the forecasts with the true values
test['unique_id'] = test['unique_id'].astype(int)
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1
ds  unique_id  y  ARCH(2)  

0  20230525  1  0.875758  1.681836 
1  20230526  1  1.304909  0.777028 
2  20230530  1  0.001660  0.968701 
…  …  …  …  … 
84  20230926  1  1.473453  NaN 
85  20230927  1  0.022931  NaN 
86  20230928  1  0.589317  NaN 
# Merge the forecasts with the true values
fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds')
plot_df[['y', "ARCH(2)"]].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Year ', fontsize=20)
ax.set_xlabel('Timestamp [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)
plt.show()
Predict method with confidence interval
To generate forecasts use the predict method.
The predict method takes two arguments: forecasts the next h
(for
horizon) and level
.

h (int):
represents the forecast h steps into the future. In this case, 12 months ahead. 
level (list of floats):
this optional parameter is used for probabilistic forecasting. Set the level (or confidence percentile) of your prediction interval. For example,level=[95]
means that the model expects the real value to be inside that interval 95% of the times.
The forecast object here is a new data frame that includes a column with the name of the model and the y hat values, as well as columns for the uncertainty intervals.
This step should take less than 1 second.
sf.predict(h=horizon)
ds  ARCH(2)  

unique_id  
1  20230525  1.681836 
1  20230526  0.777028 
1  20230529  0.677960 
…  …  … 
1  20230920  0.136752 
1  20230921  0.082173 
1  20230922  0.450958 
forecast_df = sf.predict(h=horizon, level=[80,95])
forecast_df
ds  ARCH(2)  ARCH(2)lo95  ARCH(2)lo80  ARCH(2)hi80  ARCH(2)hi95  

unique_id  
1  20230525  1.681836  0.419322  0.307963  3.055710  3.782995 
1  20230526  0.777028  3.939044  2.844559  1.290504  2.384989 
1  20230529  0.677960  3.907244  2.789475  1.433555  2.551323 
…  …  …  …  …  …  … 
1  20230920  0.136752  0.795371  0.472731  0.746235  1.068876 
1  20230921  0.082173  0.852268  0.528825  0.693172  1.016615 
1  20230922  0.450958  1.337117  1.030386  0.128471  0.435202 
We can join the forecast result with the historical data using the
pandas function pd.concat()
, and then be able to use this result for
graphing.
df_plot=pd.concat([df, forecast_df]).set_index('ds').tail(220)
df_plot
unique_id  y  ARCH(2)  ARCH(2)lo95  ARCH(2)lo80  ARCH(2)hi80  ARCH(2)hi95  

ds  
20230321  1  1.298219  NaN  NaN  NaN  NaN  NaN 
20230322  1  1.646322  NaN  NaN  NaN  NaN  NaN 
20230323  1  0.298453  NaN  NaN  NaN  NaN  NaN 
…  …  …  …  …  …  …  … 
20230920  NaN  NaN  0.136752  0.795371  0.472731  0.746235  1.068876 
20230921  NaN  NaN  0.082173  0.852268  0.528825  0.693172  1.016615 
20230922  NaN  NaN  0.450958  1.337117  1.030386  0.128471  0.435202 
def plot_forecasts(y_hist, y_true, y_pred, models):
_, ax = plt.subplots(1, 1, figsize = (20, 7))
y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
df_plot = pd.concat([y_hist, y_true]).set_index('ds').tail(12*10)
df_plot[['y'] + models].plot(ax=ax, linewidth=2 )
colors = ['green']
# Specify graph features:
ax.fill_between(df_plot.index,
df_plot['ARCH(2)lo80'],
df_plot['ARCH(2)hi80'],
alpha=.20,
color='lime',
label='ARCH(2)_level_80')
ax.fill_between(df_plot.index,
df_plot['ARCH(2)lo95'],
df_plot['ARCH(2)hi95'],
alpha=.2,
color='white',
label='ARCH(2)_level_95')
ax.set_title('', fontsize=22)
ax.set_ylabel("Return", fontsize=20)
ax.set_xlabel('MonthDays', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)
plt.show()
plot_forecasts(train, test, forecast_df, models=["ARCH(2)"])
Let’s plot the same graph using the plot function that comes in
Statsforecast
, as shown below.
Crossvalidation
In previous steps, we’ve taken our historical data to predict the future. However, to asses its accuracy we would also like to know how the model would have performed in the past. To assess the accuracy and robustness of your models on your data perform CrossValidation.
With time series data, Cross Validation is done by defining a sliding window across the historical data and predicting the period following it. This form of crossvalidation allows us to arrive at a better estimation of our model’s predictive abilities across a wider range of temporal instances while also keeping the data in the training set contiguous as is required by our models.
The following graph depicts such a Cross Validation Strategy:
Perform time series crossvalidation
Crossvalidation of time series models is considered a best practice but most implementations are very slow. The statsforecast library implements crossvalidation as a distributed operation, making the process less timeconsuming to perform. If you have big datasets you can also perform Cross Validation in a distributed cluster using Ray, Dask or Spark.
In this case, we want to evaluate the performance of each model for the
last 5 months (n_windows=5)
, forecasting every second months
(step_size=12)
. Depending on your computer, this step should take
around 1 min.
The cross_validation method from the StatsForecast class takes the following arguments.

df:
training data frame 
h (int):
represents h steps into the future that are being forecasted. In this case, 12 months ahead. 
step_size (int):
step size between each window. In other words: how often do you want to run the forecasting processes. 
n_windows(int):
number of windows used for cross validation. In other words: what number of forecasting processes in the past do you want to evaluate.
crossvalidation_df = sf.cross_validation(df=train,
h=horizon,
step_size=6,
n_windows=5)
The crossvaldation_df object is a new data frame that includes the following columns:
unique_id:
index. If you dont like working with index just run crossvalidation_df.resetindex()ds:
datestamp or temporal indexcutoff:
the last datestamp or temporal index for the n_windows.y:
true value"model":
columns with the model’s name and fitted value.
crossvalidation_df
ds  cutoff  y  ARCH(2)  

unique_id  
1  20221221  20221220  0.605272  1.889850 
1  20221222  20221220  2.492167  0.850434 
1  20221223  20221220  1.113775  0.742012 
…  …  …  …  … 
1  20230522  20230123  0.015503  0.135570 
1  20230523  20230123  1.122203  0.081367 
1  20230524  20230123  0.731860  0.446374 
Model Evaluation
We can now compute the accuracy of the forecast using an appropiate
accuracy metric. Here we’ll use the Root Mean Squared Error (RMSE). To
do this, we first need to install datasetsforecast
, a Python library
developed by Nixtla that includes a function to compute the RMSE.
!pip install datasetsforecast
from datasetsforecast.losses import rmse
The function to compute the RMSE takes two arguments:
 The actual values.
 The forecasts, in this case,
ARCH
.
rmse = rmse(crossvalidation_df['y'], crossvalidation_df["ARCH(2)"])
print("RMSE using crossvalidation: ", rmse)
RMSE using crossvalidation: 1.3816124
As you have noticed, we have used the cross validation results to perform the evaluation of our model.
Now we are going to evaluate our model with the results of the predictions, we will use different types of metrics MAE, MAPE, MASE, RMSE, SMAPE to evaluate the accuracy.
from datasetsforecast.losses import mae, mape, mase, rmse, smape
def evaluate_performace(y_hist, y_true, y_pred, model):
y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
evaluation = {}
evaluation[model] = {}
for metric in [mase, mae, mape, rmse, smape]:
metric_name = metric.__name__
if metric_name == 'mase':
evaluation[model][metric_name] = metric(y_true['y'].values,
y_true[model].values,
y_hist['y'].values, seasonality=12)
else:
evaluation[model][metric_name] = metric(y_true['y'].values, y_true[model].values)
return pd.DataFrame(evaluation).T
evaluate_performace(train, test, Y_hat, model="ARCH(2)")
mae  mape  mase  rmse  smape  

ARCH(2)  0.935485  1064.149021  NaN  1.152612  138.403076 
Acknowledgements
We would like to thank Naren Castellon for writing this tutorial.
References
 Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
 Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica: Journal of the econometric society, 9871007..
 James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.
 Nixtla Parameters.
 Pandas available frequencies.
 Rob J. Hyndman and George Athanasopoulos (2018). “Forecasting principles and practice, Time series crossvalidation”..
 Seasonal periods Rob J Hyndman.