Copying "Introduction to Statistical Analysis" in Python Vol. 21 - Chapter 7 "Introduction to Time Series Analysis" (2) 3-term Moving Average, Exponential Smoothing, and Autoregressive Model

Chapter 7 "Introduction to Time Series Analysis"

Author of the book: Dr. Sadao Ishimura

This article covers the "Introduction to Statistical Analysis" Chapter 7 "Introduction to Time Series Analysis" Python copying activity.

This is a copying series where I calmly convert the book's figures, tables, and calculations into Python code.
In Chapter 7, we will tackle time series data analysis.

This article practices moving average, exponential smoothing, and autoregressive model.
I will continue to utilize ChatGPT as well!

Now, let's open the book and set off on a journey of statistical analysis 🚀

Introduction

This blog series introduces the "joy of statistical analysis" gained through copying the Python code from the book "Introduction to Statistical Analysis" (Tokyo Tosho, referred to as "the text").

The book introduction and citation notation are posted in the linked article.

Chapter 7 Introduction to Time Series Analysis

This article covers the following sections of Chapter 7.

The data used in this article is cited directly from the data published in the text.
For data with a small number of entries, I register the data in the code, and for data with a large number of entries, I create a CSV file and load the data.

Import the libraries used in Chapter 7.

### インポート

# 数値計算
import statistics                           # Python標準ライブラリ
import numpy as np
import pandas as pd

# 統計
import scipy.stats as stats
import pymannkendall as mk                  # mann-kendall検定（トレンドの検定）

# 時系列分析
import statsmodels.api as sm
from scipy.signal import periodogram        # ピリオドグラム
import statsmodels.tsa.api as tsa           # 指数平滑法

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

Introduction

ChatGPT as a guide to the unknown in a "journey through time".

🌿Toward an analysis that gently touches the future

Time series data is engraved with events up to today.
But that is not all—
Hints about the future, such as "what will happen next," are quietly buried within it.

What I will introduce this time are methods for smoothing out the "flow" of such data,
and gently predicting what it will look like a little bit into the future.

Moving average smooths out the waves,
exponential smoothing offers a gentle glimpse into the future, and the
autoregressive model captures the "connection with itself" that the data tells.

By capturing the gentle flow, the landscape of the future gradually comes into view—
let's experience that feeling together.

➡️ Section 7.2 3-term Moving Average

According to the text, "a moving average is a method of smoothly transforming fluctuations in time series data," and "by performing a moving average, one can bring the trend of the time series data to the surface."

A moving average is something that softens (smooths) the bumps in the data.
A 3-term moving average is the "average of three time points: a certain point and the points before and after it."
In the world of stock prices, there are "5-day moving average lines," "25-day moving average lines," and so on!

So, let's borrow the data from Table 7.2.1 "Number of Patients" in the text and draw the chart for Figure 7.2.1 "Number of Patients."

### 受療者 p.274 表7.2.1

# データの登録
data1 = pd.DataFrame(
    {'受療者数': [576, 626, 754, 727, 823, 855, 766, 943, 926, 1005, 1092, 1105]},
    index=list(range(1994, 2006)))
data1.index.name = '年'

# 結果の表示
data1

[Execution Result]
12 years from 1994: This is annual data.

### データの可視化 p.274 図7.2.1
data1.plot(marker='D')
plt.xlim(data1.index.min(), data1.index.max())
plt.ylim(380, 1220);

[Execution Result]
A year-on-year increasing trend can be seen.

Let's calculate the 3-term moving average of this data.
We will use the pandas rolling() method.
Set the period to average with 'window', and set 'center=True' to place the average value at the middle time point.
You can calculate the sum of the three periods with '.sum()' and the average of the three periods with '.mean()'.

### 3項移動平均 p.275 表7.1.2

# rollingメソッドで窓の幅=3、中心に結果を格納する方法で合計を算出
data1['3項の合計'] = data1['受療者数'].rolling(window=3, center=True).sum()
# rollingメソッドで窓の幅=3、中心に結果を格納する方法で平均値を算出
data1['3項の平均値'] = data1['受療者数'].rolling(window=3, center=True).mean()
# 結果の表示
data1.round(2) 

[Execution Result]
Calculation complete for now.

Let's visualize the "smooth" time series of the 3-term moving average.

### 3項移動平均のグラフ p.275 図7.2.2

# 3項移動平均のみのグラフの描画
data1[['3項の平均値']].plot(marker='D', color='tab:orange')
plt.xlim(data1.index.min(), data1.index.max())
plt.ylim(380, 1220);

# 原系列と3項移動平均を比べるグラフの描画
data1[['受療者数', '3項の平均値']].plot(marker='D')
plt.xlim(data1.index.min(), data1.index.max())
plt.ylim(380, 1220);

[Execution Result]
The top corresponds to Figure 7.2.2 "3-term Moving Average" in the text.
The bottom is a figure showing the original time series data and the 3-term moving average side by side.

Looking at the bottom figure, the bumps in the original data are smoothed out in the 3-term moving average.
The "consistently increasing trend" becomes clearly visible.

Note that in the process of taking the average, data at the "both ends" time points will be lost.

➡️ Section 7.3 Exponential Smoothing

According to the text, "exponential smoothing can be considered a time series analysis for forecasting," and it introduces two formulas.

Exponential Smoothing Formula

① Definition of exponential smoothing
The 1-period ahead forecast value $\hat{x}(t, 1)$ at time $t$ is defined as follows, using the weight (smoothing coefficient) $\alpha$ and the observed values $x(t), x(t-1), x(t-2), \cdots$ at and before time $t$.

When the weight is $$\alpha=0.7$$, the overall weight at each time point is

• The weight at time $${t}$$ is $$\alpha=0.7$$ in total

• The weight at time $${t-1}$$ is $$\alpha(1 - \alpha)=0.7 \times 0.3=0.21$$ in total

• The weight at time $${t-2}$$ is $$\alpha(1 - \alpha)^2=0.7 \times 0.3^2=0.063$$ in total

As shown, recent time points have larger weights and thus a greater impact on the predicted value, while the weights decrease "exponentially" as we go further into the past, reducing their impact on the predicted value.

◆

(2) Recursive representation using the predicted value from one period ago

This formula is simple!
However, the prediction for time $${t-1}$$ uses the predicted value for time $${t-2}$$, and the prediction for time $${t-2}$$ uses..., requiring a "recursive"遡 back into the past.

Let's write a "one-step-ahead prediction function using exponential smoothing" using the method in (2) that uses the predicted value from one period ago.

### 指数平滑化を用いた1期先予測関数
def exponential_smoothing_pred(list_obs, alpha):

    # 1期先予測関数
    def prediction(alpha, prev_obs, prev_pred):
        return alpha * prev_obs + (1 - alpha) * prev_pred

    # numpy配列、pandasシリーズをリスト化
    list_obs = list(list_obs)
    
    # 予測値リストの初期化: 最初の2期間に1期目の観測値を設定
    list_pred = [list_obs[0]] * 2

    # 観測値の3期~最終-1期の予測を実行して予測値リストに格納
    for obs in list_obs[1:-1]:
        pred = list_pred[-1]
        list_pred.append(prediction(alpha, obs, pred))
    
    # 誤差平方和の算出
    SSE = sum([(obs - pred)**2 for obs, pred in zip(list_obs, list_pred)])
    
    # 観測値の最終期の1期先の予測を実行
    one_period_forecast = prediction(alpha, list_obs[-1], list_pred[-1])

    # 戻り値: 1期先予測値、誤差平方和、観測値に対応する予測リスト
    return one_period_forecast, SSE, list_pred

Exponential Smoothing Example p.278

Let's use the "Bail Rate" data from Table 7.3.1 in the text to practice one-step-ahead prediction using exponential smoothing.

■ Preparation for analysis
Register the data.

### 指数平滑化 保釈率 p.278 表7.3.1

# データの登録
data2 = pd.DataFrame(
    {'保釈率': [47, 55, 41, 40, 36, 44, 32, 28, 27, 26, 27, 33]},
    index=range(1, 13))
data2.index.name = '時点t'
# 結果の表示
data2

【Execution Result】
This is 12 periods of data.

Let's visualize the time series trend.
Using the ".plot" method on a pandas DataFrame results in

### データの可視化
data2.plot(marker='D', ylim=(20, 60));

【Execution Result】
Overall, it appears to be on a downward trend.
Increases can also be seen at time points $${2, 6, 11, 12}$$.

■ One-step-ahead prediction with exponential smoothing: Specifying weights
Let's perform one-step-ahead prediction using exponential smoothing.
First, we calculate with the weight $$\alpha=0.8$$ as per the text.

### 指数平滑化を用いた1期先予測 

# 関数利用

# 設定
alpha = 0.8
# 指数平滑化の実行
forecast, SSE, preds = exponential_smoothing_pred(data2['保釈率'], alpha)
# 結果の表示
print(f'第{data2.index.max()+1}期の予測値: {forecast:.6f}')
print(f'SSE: {SSE:.5f}')
pd.DataFrame(dict(予測値=preds[1:] + [forecast]), index=range(1, len(preds)+1)
            ).round(2)

【Execution Result】
We are steadily accumulating predicted values from the first period.
Otherwise, we wouldn't be able to predict the 13th period...
SSE is the sum of squared errors (squared error) between observed values and predicted values.

Let's compare the observed values and the predicted values.

### データの可視化
data2.plot(marker='D', ylim=(20, 60))
plt.plot(range(2, 14), preds[1:] + [forecast], marker='D', label='予測値')
plt.legend();

[Execution Result]
The orange predicted values look like the blue observed values shifted by one period.
Could it be that "predicting the same value as the previous period" is actually fine!?

Let's practice exponential smoothing using a Python library.
We will use statsmodels.

🖲️ statsmodels' SimpleExpSmoothing()

# statsmodels利用

# 設定
alpha = 0.8

# 指数平滑化の実行
fit = tsa.SimpleExpSmoothing(data2).fit(smoothing_level=alpha, optimized=False)

# 結果の表示
forecast = fit.forecast(1)
print(f'第{forecast.index[0]}期の予測値: {forecast.values[0]:.6f}')
print(f"α: {fit.params['smoothing_level']:.5f}")
print(f'SSE: {fit.sse:.5f}')
fit.level.rename('予測値').to_frame().round(2)

[Execution Result]
It matches the results of the custom function!

◆

■ One-period-ahead prediction with exponential smoothing: Estimating weights
Exponential smoothing in statsmodels also calculates the optimal value for the weight $$\alpha$$.
Let's try predicting right away.

# statsmodels利用 αの最適値を推定

# 指数平滑化の実行
fit = tsa.SimpleExpSmoothing(data2).fit(optimized=True)

# 結果の表示
forecast = fit.forecast(1)
print(f'第{forecast.index[0]}期の予測値: {forecast.values[0]:.6f}')
print(f"α: {fit.params['smoothing_level']:.5f}")
print(f'SSE: {fit.sse:.5f}')
fit.level.rename('予測値').to_frame().round(2)

[Execution Result]
$$\alpha \approx 0.69$$ was obtained.
The predicted value for period 13 is smaller than when $$\alpha=0.8$$.

Let's visualize it.

### データの可視化 statsmodelsのα最適値
data2.plot(marker='D', ylim=(20, 60))
plt.plot(range(2, 14), fit.level, marker='D', label='予測値')
plt.legend();

[Execution Result]
It doesn't look much different from when $$\alpha=0.8$$.

You can also search for the optimal value of the weight $$\alpha$$ using a custom function.
We search through $$\alpha = [0.0, 0.1, 0.2, \ldots, 0.9, 1.0]$$ to find the $$\alpha$$ that minimizes the sum of squared errors (SSE).

### 誤差平方和の最小となるαの探索 p.277

# 設定
alphas = np.arange(0, 1.1, 0.1)  # 探索するαの値
forecasts, SSEs = [], []         # リストの初期化

# 指数平滑化を用いた予測の実行
for alpha in alphas:
    forecast, SSE, pred_list = exponential_smoothing_pred(data2['保釈率'], alpha)
    forecasts.append(forecast)
    SSEs.append(SSE)

# 結果の表示
(pd.DataFrame({'α': alphas, '誤差平方和': SSEs, '予測値': forecasts})
 .sort_values(['誤差平方和'])
 .reset_index(drop=True)
).round(2)

[Execution Result]
The sum of squared errors was minimized when $$\alpha=0.7$$.

Understanding Check: Exponential Smoothing p.279

Using time series data for timber shipments, we perform a one-period-ahead prediction using exponential smoothing with a weight of $$\alpha=0.7$$.

I will write the code steadily.
I will borrow the data from the text.

### 指数平滑化 ブナの丸太の出荷量 p.279 表7.3.3

# データの登録
data3 = pd.DataFrame(
    {'採伐量': [2206, 2406, 2259, 2407, 2718, 2267, 2089, 1868, 1778, 1577,
               1486, 1999, 1059, 1122, 1034, 960, 938, 854, 767, 805]},
     index=range(1, 21))
data3.index.name = '時点t'

# 結果の表示
data3

[Execution Result]

Visualize the time series trend.

# 可視化
data3.plot(marker='D', ylim=(0, 2900));

[Execution Result]
A downward trend can be seen.
Increases are observed at time points $${2, 5, 12, 14}$$.

Perform a one-period-ahead prediction using exponential smoothing.
First, using the custom function.

### 指数平滑化

# 関数利用
# 設定
alpha = 0.7
# 指数平滑化の実行
forecast, SSE, preds = exponential_smoothing_pred(data3['採伐量'], alpha)
# 結果の表示
print(f'第{data3.index.max()+1}期の予測値: {forecast:.6f}')
print(f'SSE: {SSE:.5f}')

# 回答の表を作成(t=20時点まで)
data3_pred = data3.copy()
data3_pred['予測値'] = [np.nan] + preds[1:]
data3_pred.round(2)

[Execution Result]
The predicted value is $${804}$$. It feels like it reflects the downward trend.

Let's overlay the observed values and predicted values for visualization.

### データの可視化
data3.plot(marker='D', ylim=(0, 2900))
plt.plot(range(2, 22), preds[1:] + [forecast], marker='D', label='予測値')
plt.legend();

[Execution Result]
It feels like it is significantly influenced by the increases and decreases in the observed values.
It gives the impression that the observed value from the previous period is sliding into the predicted value.

We will also use statsmodels to estimate the weight $${\alpha}$$.

🖲️statsmodels' SimpleExpSmoothing()

# statsmodels利用 αの最適値を推定

# 指数平滑化の実行
fit = tsa.SimpleExpSmoothing(data3).fit()
# 結果の表示
forecast = fit.forecast(1)
print(f'第{forecast.index[0]}期の予測値: {forecast.values[0]:.6f}')
print(f"α: {fit.params['smoothing_level']:.5f}")
print(f'SSE: {fit.sse:.5f}')
fit.level.rename('予測値').to_frame().round(2)

[Execution Result]
The estimated value of $${\alpha}$$ is $${0.732}$$, and the predicted value for period 21 is $${803}$$.

Let's visualize it.

### データの可視化 statsmodelsのα最適値
data3.plot(marker='D', ylim=(0, 2900))
plt.plot(range(2, 22), fit.level, marker='D', label='予測値')
plt.legend();

[Execution Result]

➡️ Section 7.4 Autoregressive AR(1) Model

The autoregressive model is based on the idea that the observed value at a certain point in time is a weighted sum of past observed values plus white noise.
It feels like performing regression analysis using one's 'past self' as a variable, hence 'autoregressive'.

The AR(1) model is an autoregressive model that 'adds white noise to the observed value of the previous period'.
The number in parentheses is called the 'order'; in the case of AR(1), it is 'order 1'.
By the way, there are also autoregressive models of other orders, such as AR(2) and AR(3)!

I will quote the 'Definition of Time Series AR(1) Model' from the text.

When time series data is defined as $${\{\cdots, x(t-3), x(t-2), x(t-1), x(t)\}}$$,

is called an autoregressive AR(1) model.
Here, $${u(t)}$$ is white noise.
In this case, the optimal predicted value for one period ahead, $${\hat{x}(t, 1)}$$, is

.

We multiply the observed value of the previous period by the weight, which is the autoregressive coefficient $${\alpha_1}$$, and add white noise.
It is a simple formula.

Example using virtual data

Since there is no example in the text, let's create virtual data exclusive to this article and build an AR(1) model.

■ Preparation for analysis
We will create virtual data using the AR(1) model formula $${x(t)=0.72 \cdot x(t_1) + white noise}$$.
The white noise is a standard normal distribution random variable.

### 時系列データの作成

## 設定
T = 100                               # 時間の数
rng = np.random.default_rng(seed=42)  # 乱数生成器

## データの生成
# ノイズの生成：標準正規分布乱数
noise = rng.standard_normal(size=T)
# 観測値obsの作成：自己回帰 x(t) = 0.72 x(t-1) + noise(t)
obs = [noise[0]]
for t in range(1, T):
    obs.append(0.72 * obs[t-1] + noise[t])

## データのまとめ
# データフレーム化
data4 = pd.DataFrame(
    {'測定値': obs},
    index=pd.date_range(start='2016-01-01', periods=100, freq='MS'))

## データの可視化
data4.plot(figsize=(7, 4), xlabel='年月', grid=True);

[Execution Result]
This is time series data for 100 months starting from January 2016.

■ Building the model
We will build an AR(1) model.
We will use the ARIMA class from Python's statsmodels library.
Write it as follows.

### AR(1)モデルの構築

fitted_model = tsa.ARIMA(data4, order=(1, 0, 0), trend='n').fit()
fitted_model.summary()

【Execution Results】
This is the model summary.
The "coef" in the "ar.L1" row in the middle section is the estimated autoregressive coefficient, $${0.7319}$$.
It is close to the $${\alpha=0.72}$$ used for data generation.
"sigma2" is the variance of the white noise, $${0.5988}$$.
It feels a bit far from the $${\sigma^2=1^2}$$ used for data generation.

■ Future forecasting
Let's perform a forecast using the AR(1) model.
We will try a 12-month forecast. Write it as follows.

### AR(1)モデルで予測

## 予測
# 予測の実行
forecasts = fitted_model.get_forecast(steps=12)
# 予測値（平均）と信頼区間の取得
forecasts_mean = forecasts.predicted_mean
conf_ints = forecasts.conf_int()

## 描画
# 描画領域の設定
plt.figure(figsize=(7, 4))
# 観測値の描画（青色）
plt.plot(data4, label='観測値')
# 予測値（平均）の描画（赤色）
plt.plot(forecasts_mean, color='tab:red', label='予測値')
# 予測値の95%信頼区間の塗りつぶし
plt.fill_between(conf_ints.index, conf_ints.iloc[:, 0], conf_ints.iloc[:, 1], 
                 color='tomato', alpha=0.2, label='95%信頼区間')
# 観測値の最終時点に垂直点線を描画
plt.axvline(data4.index[-1], color='gray', lw=3, ls='--')
# 修飾
plt.xlabel('年月')
plt.grid()
plt.legend();

【Execution Results】
The red line is the mean of the predicted values, and the light red area is the 95% confidence interval of the predicted values.
The jagged movements of the observed values have been smoothed out.

Introduction to a book where you can learn "what comes next" in time series data

📗 Practical Data Science Series: Introduction to Time Series Analysis with Python
This is a great book for learning time series analysis in Python from basics to applications.
I recommend it to those who want to advance their learning of time series analysis.

Let's have ChatGPT wrap up the end of the article!
The "story of lines" will surely continue from here.

📘 A word from ChatGPT:

That is all for this coding session.

This article is the final installment of the series.
Thank you for staying with me for so long!

I would like to conclude with ChatGPT.

The End

Series Articles

Previous Article

Table of Contents

Blog Introduction

I am writing seven series of articles on note.
Please take a look!

1. Relaxed Statistics

This is a blog that roughly explores probability and statistics using the Statistical Test Grade 2 problem collection as a guide.
Feel free to read it as if we are just chatting. Please take a look.
It corresponds to the CBT-compatible version of the official Statistical Test Grade 2 problem collection.
Sample code for Python and EXCEL is also available.

2. Experiment! Fun Bayesian Modeling 1 & 2 with PyMC Ver. 5

I will draw and analyze the Bayesian models used in psychological research from the books 'Fun Bayesian Modeling' and 'Fun Bayesian Modeling 2' using PyMC Ver. 5.
Many Bayesian models, including those in these books, are written in R + Stan.
I strive to explore the possibilities of PyMC and make Bayesian modeling easy to practice.
Since these are familiar and easy-to-visualize themes, please try running them with PyMC and let's enjoy it together!

3. Experiment! Bayesian Modeling from Iwanami Data Science 1 with PyMC Ver. 5

I will draw and analyze the Bayesian models by four Bayesians from the book 'Experiment! Iwanami Data Science Vol. 1' using PyMC Ver. 5.
This book is a great resource for learning the basics of Bayesian programming.
I feel like I've become friends with Bayesian methods by happily running PyMC models.
Everyone, please try running them with PyMC and let's play and learn together!

4. Fun Copying: Bayesian, Python, etc.

I will blog about the results of my 'book copying activities' for Bayesian, Python, and others.
I am mainly working on translations into Python.
I hope this serves as sample code for fellow learners who are also copying code 🍀

5. Introduction to Time Series Analysis for Psychology with R and Stan, using Python and PyMC Ver. 5

I will practice the time series analysis from the book 'Introduction to Time Series Analysis for Psychology with R and Stan' using Python and PyMC Ver. 5.
This book is packed with themes on time series analysis!
I realized the depth of time series analysis.
I will enjoy learning time series analysis with my favorite language, Python.

6. Writing about Data Science-like things

I write columns on statistics, data analysis, AI, machine learning, and Python on an irregular basis.
There are many articles related to statistics and data science books.
Series on 'Statistics', 'Python', 'Mathematics and Python', and 'R' have been created.

7. Python Machine Learning Programming Practice Log

I wrote articles about my various thoughts when studying the book 'Python Machine Learning Programming: PyTorch & scikit-learn Edition'.
This book is a textbook for scikit-learn and PyTorch.
Please feel free to try it out.

Thank you very much for reading until the end.