Rによるデータサイエンス（20/21） - 「大人の教養・知識・気付き」を伸ばすブログ

　Rについて

Rによるデータサイエンス(第2版):データ解析の基礎から最新手法まで

作者:金明哲
森北出版

Amazon

をベースに学んでいく。
　今回は時系列分析（PP.269-295）を扱う*1。

前回

power-of-awareness.com

前回
21.　時系列分析
補足　スペック情報

21.　時系列分析

　時間と共に変動する現象に対して、時間の順序で測定・観測した結果を記録したデータを時系列データという。時系列データが伴う変動の振る舞いを統計的に分析するのが時系列データ分析の主要な目的で、そのような分野を時系列解析という。

21.1　時系列分析の基本概念

　時系列データは通常、ある一定の時間間隔で測定・観測したものである。観測・測定値 $y$ に関する観測回数 $n$ の時系列データは、測定時点 $t$ を添字として

$\begin{aligned} y_1,y_2,\cdots,y_{t-1},y_t,y_{t+1},\cdots,y_n \end{aligned}$

のように記述される。
　 $\mathrm{R}$ には時系列データを専門に扱うデータ型が存在する。

##################
### 時系列解析 ###
##################

library("ggplot2")
library("ggfortify")

# データの読込
data(lh)
data(UKgas)

class(lh) # データ型
print(lh) # データ内容
print(UKgas)

# 開始時点と終了時点
start(UKgas)
end(UKgas)

# 特定期間を抜き出す
window(UKgas, c(1975,2),c(1979,3)) # 1975年2Qから1979年3Qまで


## 時系列データの図示
g_ts <- autoplot(UKgas)
g_ts <- g_ts + theme_classic()
g_ts <- g_ts + labs(x="年(四半期)", y = "天然ガス消費量[百万サーム]")
g_ts <- g_ts + theme(plot.title = element_text(hjust = 0.5, size = 14),
                     axis.title = element_text(size = 14),
                     axis.text = element_text(size = 14),
                     legend.title = element_text(size = 14),
                     legend.text = element_text(size = 14),
                     legend.position = "none")
g_ts <- g_ts + scale_y_continuous(breaks = seq(0,1200,100),labels = seq(0,1200,100))
g_ts <- g_ts + geom_hline(yintercept = 0, linetype = 1)
g_ts <- g_ts + geom_hline(yintercept = seq(0,1200,100), linetype = 3)

plot(g_ts)

21.1.1　データオブジェクトの作成

### tsオブジェクトの作成

tmp <- ts(1:120,c(1995,6), frequency = 12)
class(tmp)

print(tmp)

21.1.2　ラグ処理

## ラグ処理
library("datasets")

ldeaths
lag(ldeaths,k = 1)

plot(diff(UKgas,lag = 4))

# 以下は等価
diff(UKgas, lag = 2)
UKgas - lag(UKgas, k = -2)

21.2　自己共分散と自己相関

　時系列 $\{y_1,\cdots,y_n\}$ に対して各種統計量を定義できる。標本平均、自己共分散関数および自己相関関数はそれぞれ

$\begin{aligned} \hat{\mu}&=\displaystyle{\frac{1}{n}\sum_{t=1}^{n}y_t},\\ \hat{C}_k&=\displaystyle{\frac{1}{n}\sum_{t=k+1}^{n}(y_t-\hat{\mu})(y_{t-k}-\hat{\mu}) },\\ \hat{\rho}_k&=\displaystyle{\frac{\hat{C}_k}{\hat{C}_0}} \end{aligned}$

で定義される。
　パッケージ $\mathrm{stats}$ には自己相関および自己共分散を計算する関数 $\mathrm{acf}$ がある。

### 自己共分散と自己相関
library("stats")

ls_acf <- acf(UKgas)

df_acf <- data.frame(lag = ls_acf$lag,
                     acf = ls_acf$acf)

# 
g1 <- ggplot(df_acf,aes(x = lag, y = acf),type = "correlation")
g1 <- g1 + geom_bar(stat = "identity")
g1 <- g1 + theme_classic()
g1 <- g1 + labs(x="ラグ", y = "標本自己相関関数")
g1 <- g1 + theme(plot.title = element_text(hjust = 0.5, size = 14),
                     axis.title = element_text(size = 14),
                     axis.text = element_text(size = 14),
                     legend.title = element_text(size = 14),
                     legend.text = element_text(size = 14),
                     legend.position = "none")
g1 <- g1 + scale_y_continuous(breaks = seq(0,1,0.25),labels = seq(0,1,0.25))
g1 <- g1 + geom_hline(yintercept = 0, linetype = 1)
g1 <- g1 + geom_hline(yintercept = seq(0,1,0.25), linetype = 3)

plot(g1)

21.3　スペクトル分析

　時系列データは、トレンド、周期成分およびノイズなど、複数の成分が混合しているものと考えられる。時系列データに隠された周期性を解析する方法としてスペクトル解析がある。
　時系列における自己共分散 $C_k$ の $\mathrm{Fourier}$ 変換が可能であるとき、 $-\displaystyle{\frac{1}{2}}\leq f\leq\displaystyle{\frac{1}{2}}$ 上で定義される関数

$\begin{aligned} p(f)=\displaystyle{\sum_{k=-\infty}^{\infty}C_k\exp\left(-2\pi i kf \right)}=C_0+2\displaystyle{\sum_{k=1}^{\infty}C_k\cos(2\pi kf)} \end{aligned}$

をパワースペクトル密度関数と呼ぶ。標本 $y_1,\cdots,y_n$ の自己共分散 $\hat{C}_k$ を用いて定義したもの

$\begin{aligned} p_j=\hat{C}_0+2\displaystyle{\sum_{k=1}^{n-1}\hat{C}_k\cos(2\pi kf)} \end{aligned}$

をピリオドグラムという。ただし周波数は $f_j=\displaystyle{\frac{j}{n}},j=0,1,\cdots,\displaystyle{\frac{n}{2}}$ である。

library("ggplot2")

data(UKgas)
vc_frq <- spec.pgram(UKgas)$freq
vc_spc <- spec.pgram(UKgas)$spec

df_spc <- data.frame(frq = vc_frq,
                     var = "spectrum",
                     spc = log10(vc_spc))

g2 <- ggplot(df_spc,aes(x = frq, y = spc, color = var, fill = var))
g2 <- g2 + geom_line() + geom_point()
g2 <- g2 + theme_classic()
g2 <- g2 + labs(x = "周波数", y = "スペクトル")
g2 <- g2 + theme(plot.title = element_text(hjust = 0.5, size = 14),
                 axis.title = element_text(size = 14),
                 axis.text = element_text(size = 14),
                 legend.title = element_text(size = 14),
                 legend.text = element_text(size = 14),
                 legend.position = "none")
g2 <- g2 + scale_x_continuous(breaks = seq(0,vc_frq[length(vc_frq)],0.5),
                              labels = seq(0,vc_frq[length(vc_frq)],0.5))
g2 <- g2 + scale_y_continuous(breaks = seq(floor(min(log10(vc_spc))),ceiling(max(log10(vc_spc))),1),
                              labels = 10^seq(floor(min(log10(vc_spc))),ceiling(max(log10(vc_spc))),1))
g2 <- g2 + geom_hline(yintercept = 0, linetype = 1)
g2 <- g2 + geom_hline(yintercept = seq(floor(min(log10(vc_spc))),ceiling(max(log10(vc_spc))),1), linetype = 3)

plot(g2)

21.4　ランダムウォークと単位根

　時系列解析も、回帰分析と同じく、収集したデータからモデルを構築し、将来の予測やシステム制御等を行うのを目的とする。ある時系列データが、

$\begin{aligned} y_{t+1}=ay_{t}+\varepsilon_{t+1} \end{aligned}$

で表せ、更に単位根、すなわと $|a|=1$ のとき、この $y_t$ はランダムウォークであるという。ランダムウォークは非定常であるから、時系列データを分析するには、まずそのデータがランダムウォークであるか否かを調べることが重要である。
　時系列データが単位根なのかを検定することを単位根検定という。すなわち単位根検定は

帰無仮説 $H_0$ ：単位根が存在する(すなわち $|a|=1$ )
対立仮説 $H_1$ ：単位根が存在しない(すなわち $|a|\neq1$ )

を検定する仮説検定である。これを具体的に解くための手法はいくつか存在する。
　パッケージ $\mathrm{stats}$ には $\mathrm{Phillips}$ - $\mathrm{Perron}$ 検定が実装されている。

##################
### 単位根検定 ###
##################
library("tseries")

df_UKgas_log_diff <- log(UKgas) - log(lag(UKgas,k = 4))

# 単位根か否か
PP.test(df_UKgas_log_diff) # 帰無仮説は棄却できない=単位根とは言えない
adf.test(df_UKgas_log_diff) # 帰無仮説は棄却できない=単位根とは言えない

21.5　ARモデル

　時点 $t-p,p\gt0$ から時点 $t$ までの各データの関係式

$\begin{aligned} y_t=\displaystyle{\sum_{i=11}^{p}a_iy_{t-i}}+\varepsilon_t \end{aligned}$

を自己回帰モデル( $\mathrm{AR}$ モデル)という。
　 $\mathrm{R}$ では $\mathrm{stats}$ の $\mathrm{ar}$ ()関数を用いればよい。

# ARモデル
ls_ar_UKgas <- ar(df_UKgas_log_diff) # 特に指定しないとAICを用いて次数を自動的に決定

21.5.1　モデルの診断

　モデルの診断として、残差分析には、 $\mathrm{Box}$ - $\mathrm{Pierce}$ 検定などにより各残差間の独立性を検定できる。また $\mathrm{Jarque}$ - $\mathrm{Bera}$ 検定により、正規性の有無も検定できる。

Box.test(x = ls_ar_UKgas$resid,type = "Ljung") # Ljung-Box検定。デフォルトではBox-Pierce検定。帰無仮説を棄却できそう。
jarque.bera.test(x = ls_ar_UKgas$resid[!is.na(ls_ar_UKgas$resid)]) # Jarque-Bera検定。正規性も否定はできなさそう

21.5.2　予測

predict(ls_ar_UKgas,n.ahead = 12)

21.6　ARMA/ARIMAモデル

　 $\mathrm{AR}$ モデルに誤差の移動平均を加えたモデル

$\begin{aligned} \displaystyle{\sum_{i=1}^{p} a_iy_{t-i}}+\varepsilon_t+\displaystyle{\sum_{j=1}^{q}b_j\varepsilon_{t-j}} \end{aligned}$

を $\mathrm{ARMA}$ モデルという。
　 $y_t$ の $d$ 階差分 $\Delta^{d}y_t$ の $\mathrm{ARMA}$ モデルを自己回帰和分移動平均( $\mathrm{ARIMA}$ )モデルという。

###################
### ARIMAモデル ###
###################

ls_arima <- arima(x = log(UKgas),order = c(4,4,1))

tsdiag(ls_arima,gof.lag = 12) # tsdiag:残差分析のための関数

21.7　その他のモデル

21.7.1　ARFIMAモデル

　差分の階数 $d$ を実数に拡張したモデルを自己回帰実数和分移動平均( $\mathrm{ARFIMA}$ )モデルという。 $\mathrm{fracdiff}$ パッケージを用いると、適切な差分階数を推定した上で、モデルの係数を決めて同モデルを構築できる。

####################
### ARFIMAモデル ###
####################
library("fracdiff")

data(AirPassengers)

d_estimate <- fdGPH(AirPassengers)$d # 差分パラメータの推定

ls_arfima <- fracdiff(x = AirPassengers, nar = 3, dtol = d_estimate, nma = 1) # モデルパラメータの推定

21.7.2　GARCHモデル

　時系列データが、条件付き平均 $g_t$ および条件付き分散 $h_t$ を持つ正規分布 $N(g_t,h_t)$ に従うとし、また

$\begin{aligned} h_t=\omega+\displaystyle{\sum_{i=1}^{q}\alpha_i\varepsilon_{t-i}^2} \end{aligned}$

と書けるとするモデルを $\mathrm{ARCH}$ モデルという。これを更に拡張した

$\begin{aligned} h_t=\omega+\displaystyle{\sum_{i=1}^{q}\alpha_i\varepsilon_{t-i}^2}+\displaystyle{\sum_{j=1}^{r}\beta_i h_{t-j}} \end{aligned}$

を $\mathrm{GARCH}$ モデルという。これを更に拡張したモデルも他に提案されている。
　 $\mathrm{GARCH}$ モデルの構築には $\mathrm{tseries}$ の $\mathrm{garch}$ ()関数などを用いればよい。

21.8　成分の分解

　時系列データのうち周期性を持つものについて、

$\begin{aligned} 観測値=トレンド+周期変動+残差 \end{aligned}$

といった形で分解することで周期成分を抽出するモデルを考えることもできる。
　 $\mathrm{stats}$ パッケージの $\mathrm{stl}$ ()関数でこのモデルを考えることができる。

######################
### 周期成分の抽出 ###
#####################
df_UKgas_stl <- stl(UKgas,s.window = "per")

plot(df_UKgas_stl)

21.9　多変量時系列解析

　多変量時系列データに関するモデルを構築することもでき、その中でも $\mathrm{VAR}$ モデルが最もよく用いられる。

#################
### VARモデル ###
#################

library("tseries")
data("USeconomic")

USe_d <- diff(USeconomic[,1:2])
ts.plot(USe_d)

acf(USe_d)

USe_d_var <- ar(x = USe_d, order.max = 2, aic = F)

str(USe_d_var)
summary(USe_d_var)

plot(USe_d_var$resid)

21.10　カオス時系列

　不規則に変動する時系列データを非線形的に扱う手法として、カオス理論に基づく方法がある。そのためには $\mathrm{tseriesChaos}$ パッケージを用いる。

補足　スペック情報

エディション	Windows 10 Home
バージョン	20H2
プロセッサ	Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50 GHz
実装 RAM	8.00 GB
システムの種類	64 ビットオペレーティングシステム、x64 ベースプロセッサ
R　バージョン	4.1.3 (March, 2022)
RStudio　バージョン	2022.02.2 Build 485

*1:時系列解析は別途しっかりと扱うつもりであるため、ここでは簡単に扱う。

前回

21. 時系列分析

21.1 時系列分析の基本概念

21.1.1 データオブジェクトの作成

21.1.2 ラグ処理

21.2 自己共分散と自己相関

21.3 スペクトル分析

21.4 ランダムウォークと単位根

21.5 ARモデル

21.5.1 モデルの診断

21.5.2 予測

21.6 ARMA/ARIMAモデル

21.7 その他のモデル

21.7.1 ARFIMAモデル

21.7.2 GARCHモデル

21.8 成分の分解

21.9 多変量時系列解析

21.10 カオス時系列

補足 スペック情報