時系列分析の基礎

慶應大学 長倉大輔 様「 時系列分析 1 」を参照させていただきました。

一定の時間間隔で観測されたデータ \(y_t\)を時系列データと言います。
データを確率変数と見れば、期待値、分散が考えられます。
　　期待値：\(E(y_t)=\mu_t\)　分散：\(V(y_t)=\sigma_t^2\)

データが統計的に分析可能であるために、これらが時間tに依存しないという前提を設けます。
　　期待値：\(E(y_t)=\mu\)　分散：\(V(y_t)=\sigma^2\)

平均・分散の推定

時点 \(t=1,2,\cdots,T\) までのデータの集合 \(\left(y_1,y_2,\cdots,y_T\right)\) に対し、
　　　(標本平均) 　\(\displaystyle \overline{y}=\dfrac{1}{T}\sum ^{T}_{t=1}y_{t}\)
定常性という仮定のもとで、\(\bar{y}\) は\(E(y_t)=\mu\) の推定値となります。
時間tによって \(y_t\)の期待値が変わらないから、T個のデータから計算した平均が他の tでも当てはまると考えるのです。

　　　(分散)　　 \(V(y)=\dfrac{1}{T}\displaystyle \sum ^{T}_{t=1}\left( y_{t}-\overline{y}\right)^2\)　　も \(V(y_t)=\sigma^2\) の推定値となります。
さらに、
\(y_t\)からkだけ時間をずらした系列\(y_{t-k}\)も定常性の仮定を満たし、
　　　期待値：\(E(y_{t-k})=\mu\)　分散：\(V(y_{t-k})=\sigma^2\)

\(y_t,~y_{t-k}\)の共分散と相関係数の推定値を次のように求めることができます。

(標本自己共分散)　\(\widehat{\gamma }_{k}=\dfrac{1}{T}\sum ^{T}_{t=k+1}\left( y_{t}-\overline{y}\right) \times \left( y_{t-k}-\overline{y}\right)\)　

　(標本自己相関) 　\(\widehat{\rho}_{k}=\dfrac{\widehat{\gamma}_{k}}{\widehat{\gamma}_{0}}\)

実際の計算はRの時系列パッケージで行うことができます。

分析に必要なパッケージを読み込みます。

library(TSstudio)
library(MLmetrics)
## 
##  次のパッケージを付け加えます: 'MLmetrics'
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      Recall
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TTR)
library(quantmod)
##  要求されたパッケージ xts をロード中です
##  要求されたパッケージ zoo をロード中です
## 
##  次のパッケージを付け加えます: 'zoo'
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      as.Date, as.Date.numeric

時系列データ

時系列データとはどのようなものでしょうか。
当サイト情報Ⅰ_試作問題_参考 に時系列データの例があります。
時系列データは次の成分に分解されます。（時系列分解）

トレンド (Trend) データの長期的変化・傾向
季節変動 (Seasonal) 周期変動
誤差変動 (Random) 誤差による変化

これらの基礎成分から時系列データを合成してみます。

トレンド

0.2ずつの増加をトレンドとします。

trend <- seq(0, 8-0.2, 0.2)
str(trend)
##  num [1:40] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 ...
plot(trend,type='l')

季節成分

簡単のため8ヶ月周期の変化 ５サイクルを季節変動とします。

kisetu <- rep(c(0,2,4,2,0,-2,-4,-2), 5)
str(kisetu)
##  num [1:40] 0 2 4 2 0 -2 -4 -2 0 2 ...
plot(kisetu,type='l')

誤差成分

誤差は、平均０、標準偏差１の正規分布に従うものとします。
誤差変動は標準偏差の２倍、即ち \(\pm1\)の範囲に収まります。

gosa <- rnorm(40,0,1)
str(gosa)
##  num [1:40] 0.21 -1.122 1.235 0.626 0.608 ...
plot(gosa,type='l')

合成した時系列データ
トレンド＋季節成分＋誤差成分 を時系列データとして図示します。

x <- trend+kisetu+gosa
plot(x,type='l')

データ形式を時系列に変換
xを関数tsによって時系列データiceに変換し、パケージの関数ts_decompos で時系列分解します。

ice<-ts(x,start=c(2016,1),frequency=8)
ts_decompose(ice)

自己相関関数(ACF)

自己相関とは、\(y_t\)と \(y_t\)からkだけ時間をずらした系列\(y_{t-k}\)との相関係数です。
ACF(autocorrelation function)のグラフはラグkを横軸にとった自己相関係数のグラフです。
季節変動（12ヶ月周期）ならば、k=12で１回目のピーク、k=24で２回目のピークが観測されます。
iceの例は8ヶ月周期に作られているので、k=8ごとにピークを観測します。

str(ice)
##  Time-Series [1:40] from 2016 to 2021: 0.21 1.08 5.64 3.23 1.41 ...
acf(x,plot=T)

パッケージTS の関数ts_corを用いると、季節変動の周期（Seasonal Lag)が分析により表示され、PACFは片自己相関係数(partial autocorrelation function)で、間の観測値の影響を除去した相関係数です。

ts_cor(ice, lag.max = 60)

季節変動の大きさを変えると時系列分解はどうなるでしょうか？ 参考までに、少し変えたデータで時系列の合成と分解を行ってみます。

tre <- seq(0, 32-0.2, 0.2)
kisetu <- rep(c(0,4,8,4,0,-4,-8,-4), 20)
gosa <- rnorm(160,0,1)
x <- tre+kisetu+gosa
plot(x,type='l')

ice<-ts(x,start=c(2016,1),frequency=8)
ts_decompose(ice)
ts_cor(ice, lag.max = 60)