Hyndman and Athanasopoulos (2013) outline the general approach for fitting an ARIMA model to a set of time series data.

Plot the data. Identify any unusual observations.

If necessary, transform the data (e.g. log transformations, power transformations, etc.) to stabilize the variance.

If the data is non-stationary: take first differences of the data until the data are stationary.

Examine the ACF/PACF: Which AR($p$), MA($q$), ARMA($p, q$) or ARIMA($p, d, q$) model is appropriate?

Choose a model and then use the $AICc$ as a criterion to search for a better model.

Check the residuals from the chosen model by plotting the ACF of the residuals, and doing a

*Portmanteau test*of the residuals. If they do not look like white noise, try a modified model.Once the residuals look like white noise, calculate forecasts.

In [2]:

```
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import random
```

Examplarily, we will fit a model to the *global annual temperature anomalies* (`t_global`

) provided by the Berkeley Earth Surface Temperature Study. Temperatures are given in Celsius and are reported as anomalies relative to the period January 1850 to December 2000 average. Revisit the section on *data sets used* to remind yourself how we downloaded and extracted the data set of interest.

In [3]:

```
t_global = pd.read_json(
"http://userpage.fu-berlin.de/soga/soga-py/300/307000_time_series/t_global.json"
)
t_global["Date"] = pd.to_datetime(t_global["Date"], format="%Y-%m-%d", errors="coerce")
t_global = t_global.set_index("Date")["Monthly Anomaly_global"]
t_global
```

Out[3]:

Date 1750-01-01 -0.993 1750-02-01 -1.679 1750-03-01 -0.192 1750-04-01 -0.531 1750-05-01 -1.881 ... 2022-05-01 1.023 2022-06-01 1.315 2022-07-01 1.289 2022-08-01 1.231 2022-09-01 1.090 Name: Monthly Anomaly_global, Length: 3273, dtype: float64

We see that the data set consist of monthly mean temperatures.

Exercise: Aggregate the data (

`t_global`

) to annual means using the`groupby()`

function.

In [4]:

```
## Your code here ...
```

In [5]:

```
temp_global_year = t_global.groupby(t_global.index.to_period("Y")).agg("mean")
```

Then, we split the data set into two parts. One part, the *training set*, is the period 1850 to 2000. Based on that data set we learn the model parameters. The other part, the *test set*, is the period 2001 to 2021. This part is used to validate our model, as we compare the model forecast with the observations.

Exercise: Subset the time series (

`temp_global_year`

) to the period 1850 to 2000 and to the period 2001 to 2016.

In [6]:

```
## Your code here...
```

In [7]:

```
temp_global_training = temp_global_year["1850-01-01":"2000-01-01"]
temp_global_test = temp_global_year["2000-01-01":]
```

Let us plot the data and identify any unusual observations.

In [8]:

```
plt.figure(figsize=(18, 6))
plt.title("Earth Surface Temperature Anomalies", fontsize=14)
temp_global_training.plot(label="training set 1850-2000", fontsize=14)
temp_global_test.plot(label="test set 2001-2016", fontsize=14)
plt.legend()
plt.show()
```

It seems that there are no outliers or unusual observation in the data set.

With respect to the plot above it does not seem that the variation increases or decreases with the level of the time series. Thus, probably no transformation is required. However, just as a sanity check we apply a *Box-Cox transformation* using the `scipy.stats.boxcox`

function in Python and plot the transformed data against the original data.

The Box-Cox transformation is defined as follows:

$$ w_t = \begin{cases} log(y_t) & \text{if } \lambda=0 \\ (y_t^{\lambda}-1)/\lambda & \text{otherwise} \end{cases} $$The Box-Cox transformation depends on the parameter $\lambda$. So if $\lambda =0$ the natural logarithm is used, but if $\lambda \ne0$, a power transformation is used followed by some simple scaling. We can make use of the `scipy.stats.boxcox`

function, which chooses the value of $\lambda$, which minimizes the coefficient of variation and uses that value within the transformation.

In [9]:

```
t_global
```

Out[9]:

In [10]:

```
from scipy.stats import boxcox
boxcox_transformed_data, boxcox_lamba = boxcox(temp_global_training + 10)
boxcox_transformed_data = pd.Series(
boxcox_transformed_data, index=temp_global_training.index
)
```

We use the given $\lambda$ and plot the transformed time series against the original time series.

In [11]:

```
fig, ax = plt.subplots(2, 1, figsize=(16, 8))
temp_global_training.plot(ax=ax[0], color="black", fontsize=14)
ax[0].set_title("Original time series", fontsize=14)
boxcox_transformed_data.plot(
ax=ax[1],
color="grey",
)
ax[1].set_title("Box-Cox transformed time series", fontsize=14)
ax[0].grid()
ax[1].grid()
plt.tight_layout()
plt.show()
```

The Box-Cox transformed time series (gray line) does not differ significantly from the original time series (black line). Hence, we continue our analysis with the original time series data.

Fitting an ARIMA model requires the series to be stationary, which means that the *mean*, *variance*, and *autocovariance* of the series are time invariant. Based on our domain knowledge and the visualization of the time series data we are quite sure that the temperature time series is not stationary.

As a proof of concept we apply the **Kwiatkowski-Phillips-Schmidt-Shin (KPSS)** test. Recall, the null hypothesis in a KPSS test is that a time series is stationary. We use the function we introduced in a previous section.

$H_0:\text{The time series is trend-stationary.}$

$H_A:\text{The time series is not trend-stationary.}$

In [12]:

```
# KPSS test
from statsmodels.tsa.stattools import kpss
def kpss_test(series, **kw):
statistic, p_value, n_lags, critical_values = kpss(series, **kw)
# Format Output
print(f"KPSS Statistic: {statistic}")
print(f"p-value: {p_value}")
print(f"num lags: {n_lags}")
print("Critial Values:")
for key, value in critical_values.items():
print(f" {key} : {value}")
print(f'Result: The series is {"not " if p_value < 0.05 else ""}stationary')
```

In [13]:

```
kpss_test(temp_global_training)
```

KPSS Statistic: 1.613241516211716 p-value: 0.01 num lags: 8 Critial Values: 10% : 0.347 5% : 0.463 2.5% : 0.574 1% : 0.739 Result: The series is not stationary

C:\Users\mceck\miniconda3\envs\rasterdata\lib\site-packages\statsmodels\tsa\stattools.py:2018: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is smaller than the p-value returned. warnings.warn(

The *p*-value of the KPSS test on the original data set (`temp_global_training`

) is $p<0.01$. Hence we reject $H_0$ in favor of $H_A$; the original time series is **not** trend-stationary.

A non-stationary series can be corrected by a simple transformation such as differencing. Recall that the difference $(\Delta y_t)$ is calculated by subtracting one period's values from the previous period's values. Let us apply the `diff()`

function and compute the KPSS test for the differenced time series.

Exercise: Difference the

`temp_global_training`

time series and apply a KPSS test on the differenced time series.

In [14]:

```
## Your code here...
```

In [15]:

```
temp_global_training_diff1 = temp_global_training.diff()
```

In [16]:

```
kpss_test(temp_global_training_diff1.dropna()) ## ignore NaN for kpss
```

KPSS Statistic: 0.0976339172178381 p-value: 0.1 num lags: 24 Critial Values: 10% : 0.347 5% : 0.463 2.5% : 0.574 1% : 0.739 Result: The series is stationary

C:\Users\mceck\miniconda3\envs\rasterdata\lib\site-packages\statsmodels\tsa\stattools.py:2022: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. warnings.warn(

The *p*-value of the KPSS test on the differenced data set (`temp_global_training_diff1`

) is $p>0.1$. Hence we would not reject $H_0$ in favor of $H_A$; the differenced time series is trend-stationary.

Let us plot the differenced data set.

In [17]:

```
plt.figure(figsize=(18, 6))
plt.title("Differenced data set: temp_global_training_diff1`")
temp_global_training_diff1.plot(color="black")
plt.grid()
plt.show()
```

`plot_acf()`

function from `statsmodels.graphics.tsaplots`

.

In [18]:

```
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf
plt.figure(figsize=(18, 6))
plot_acf(temp_global_training_diff1.dropna())
plt.title("ACF for Differenced Series")
plt.show()
```

<Figure size 1800x600 with 0 Axes>