第12講時系列解析

時系列の基本モデル

Published

December 26, 2025

準備

以下で利用する共通パッケージを読み込む．

library(conflicted) # 関数名の衝突を警告
  conflicts_prefer(   # 優先的に使う関数を指定
    dplyr::filter(),
    dplyr::lag(),
    dplyr::select(),
    scales::discard(),
    recipes::fixed(),
    yardstick::spec(),
    recipes::step(),
    )
library(tidyverse)  
library(ggfortify)
library(ggrepel)
library(GGally)
library(gt)         
library(gtsummary)  
library(fable)
library(tsibble)
library(feasts)
library(patchwork)
#' macOSのための日本語表示の設定
if(Sys.info()["sysname"] == "Darwin") { # macOSか調べる
    jp_font <- "HiraMaruProN-W4"
    theme_update(text = element_text(family = jp_font))
    update_geom_defaults("text", list(family = theme_get()$text$family))
    update_geom_defaults("label", list(family = theme_get()$text$family))
    update_geom_defaults("text_repel", list(family = theme_get()$text$family))
    update_geom_defaults("label_repel", list(family = theme_get()$text$family))
} else {jp_font <- NULL}

時系列解析に用いる関数

時系列解析のためのパッケージは多岐に渡るが，本講義では tidyverse の考え方に沿って時系列解析用に開発されているパッケージ群である tidyverts を主に取り上げる．以下に本講義で扱う範囲内の base R および tidyverts の基本的な関数を説明する．

stats : base R の基本的な統計に関するパッケージ
- 関数 ts(), acf() など
- 標準でインストールされている
fable/tsibble/feasts : 時系列のためのパッケージ
- https://tidyverts.org (開発元)
- forecast パッケージの tidyverse 版
- 関数 ARIMA(), ACF(), autoplot() など

1次元の時系列はベクトル，多次元の時系列はデータフレームとして表すことができるが，時刻や周期といった時系列特有の情報を扱いやすくするために時系列特有のクラスがいくつか定義されている．最も単純な時系列の作成は以下の関数 stats::ts() を用いる．

時系列の作成

ts(data = NA, start = 1, end = numeric(), frequency = 1,
   deltat = 1, ts.eps = getOption("ts.eps"),
   class = if(nseries > 1) c("mts", "ts", "matrix", "array") else "ts",
   names = )
#' data: ベクトル，または行列(データフレーム)
#' start: 開始時刻
#' end: 終了時刻
#' frequency: 単位時間あたりの観測回数
#' 詳細は '?stats::ts' を参照

典型的な使い方は以下のようになる．

x <- rnorm(24) # 正規分布のホワイトノイズ
ts(data = x) # t=1,2,... を添字とする単純な時系列
ts(data = x, start = c(2020,1), frequency =12) # 2020年1月からの月ごと
ts(data = x, start = c(2020,3), frequency =4) # 四半期ごと

ts クラスのオブジェクトは通常その時間情報を利用して処理が行われるため，関数によっては扱いがベクトルと異なる場合があるので注意する必要がある．

tidyverts では関数 tsibble::tsibble() を用いて時系列版のtibble クラスを作成する．

時系列の作成

tsibble(..., key = NULL, index, regular = TRUE, .drop = TRUE)
#' ...: データ
#' key: indexの補助情報(同じ時間の異なるデータを表す)
#' index: 時間情報を表す列を設定
#' 詳細は '?tsibble::tsibble' を参照

ts クラスや tibble クラスを tsibble クラスに変換するには以下を用いる．

as_tsibble(x, key = NULL, index,
           regular = TRUE, validate = TRUE, .drop = TRUE, ...)
#' x: データ(時系列オブジェクトやデータフレーム)
#' 詳細は '?tsibble::as_tsibble' を参照

以下は使い方の例である．

tsibble(date = as_date("2024-01-01") + 0:9,
        value = rnorm(10))
tibble(year = 2001:2020,
       value = rnorm(20)) |>
    as_tsibble(index = year) # yearを時間情報に指定
AirPassengers |> as_tsibble() # 時系列オブジェクトの変換

時系列の視覚化には関数 fabletools::autoplot() を利用すればよい．

時系列の描画

autoplot(object, .vars = NULL, ...)
#' 詳細は '?fabletools::autoplot.tbl_ts()'を参照

以下のように1次元でも多次元でも使うことができる．

#' 単一時系列の描画
x <- rnorm(240) # 正規分布のホワイトノイズ
ts(x, start = c(2000,1), frequency = 12) |> # 2000年1月から毎月のデータ
    as_tsibble() |> # tsibbleクラスに変換
    autoplot(value) # value (as_tsibbleによる時系列データの列名の規定値) を描画
#' 複数の系列を表示する場合
y <- rt(240, df=4)  # t分布のホワイトノイズ
z <- ts(tibble(x,y), start = c(2000,1), frequency = 12)
z |> as_tsibble() |>
    autoplot(value) # 同一のグラフで色を変えて描画
z |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) + 
    facet_grid(key ~ .) # 別にする場合は facet を指定すればよい

基本的な時系列モデル

問題

指定された確率過程を生成して図示しなさい．

平均0，分散4の正規分布に従うホワイトノイズ．
上記のホワイトノイズに初期値-1で単位時刻あたり1/20で増加するトレンドを持つ確率過程．
上記のホワイトノイズから生成されるランダムウォーク．

解答欄

解答例

以下で作成する時系列に共通な変数を定義しておく．

Tmax <- 200 # 時系列の長さ t=1,..,Tmax
K <- 5 # 生成する時系列の数
#' set.seed(123) # 必要なら乱数のシード値を指定する

ホワイトノイズの生成と図示を行う．

x <- ts(rnorm(Tmax, sd=2)) # 分散4=標準偏差2
#' x <- ts(rt(n, df=4)) # 正規分布ではなく，例えば自由度4のt分布の場合
x |>
    as_tsibble() |>
    autoplot(value)

色や線種を指定することもできる．

x |>
    as_tsibble() |>
    autoplot(value,
             colour = "blue",
             linetype = "dotted")

ggplot オブジェクトとして修飾することも可能である．

x |>
    as_tsibble() |>
    autoplot(value, colour = "blue") +
    labs(x = "time", y = "observation")

トレンドのあるホワイトノイズの生成と図示を行う．

x <- ts(rnorm(Tmax, sd=2) -1 + 0.05*(1:Tmax))
x |> as_tsibble() |> autoplot(value)

ランダムウォークは定義に則って再帰的に計算すればよい．

x <- ts(rnorm(Tmax, sd=2)) # はじめは epsilon が入っている
for(t in 2:Tmax) {
    x[t] <- x[t-1] + x[t] # 順に足し合わせていく
}
x |> as_tsibble() |> autoplot(value)

同じ演算をする関数が用意されている．

x <- ts(cumsum(rnorm(Tmax, sd=2))) # 逐次的に加算を行う関数
x |> as_tsibble() |> autoplot(value)

これ以外にも書き方はいろいろあるので考えてみて下さい．

複数の系列を表示する場合には以下のようにすればよい．

ホワイトノイズの場合でいくつかの図示の方法を例示する．

z <- ts(replicate(K, # K回以下の関数を実行する
                  rnorm(Tmax, sd = 2)))
z |> as_tsibble() |>
    autoplot(value) # 同じグラフに重ね描きする

z |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) # 凡例を削除

z |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    facet_grid(key ~ .) # 各系列を個別に表示する

z |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    labs(title = expression(X[t] == epsilon[t]), # 数式で表示
         x = "Time", y = "Observations")

トレンドのあるホワイトノイズの一例は以下のようになる．

z <- ts(replicate(K,
                  rnorm(Tmax, sd = 2) -1 + 0.05*(1:Tmax)))
z |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    labs(title = expression(X[t] == -1 + 0.05 * t + epsilon[t]),
         x = "Time", y = "Observations")

ランダムウォークの一例は以下のようになる．

z <- ts(replicate(K,
                  cumsum(rnorm(Tmax, sd=2))))
z |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    labs(title = expression(X[t] == X[t-1] + epsilon[t]),
         x = "Time", y = "Observations")

より一般の時系列モデル

問題

平均0，分散1のホワイトノイズを用いて，以下の指定された確率過程を生成し，図示しなさい．

係数 a_{1}=0.67,a_{2}=0.26 を持つAR(2)過程．
係数 b_{1}=0.44,b_{2}=0.08 を持つMA(2)過程．
係数 a_{1}=0.8,a_{2}=-0.64,b_{1}=-0.5 を持つARMA(2,1)過程．

解答欄

解答例

前の練習問題と同じように設定する．

Tmax <- 200 # 時系列の長さ t=1,..,Tmax
K <- 5 # 生成する時系列の数

AR(2)モデルのシミュレーションを行う．

a <- c(0.67, 0.26) # ARの係数
epsilon <- rnorm(Tmax) # epsilonを生成
x <- double(Tmax)      # 変数を用意
x[1:2] <- epsilon[1:2] # 初期値は(epsilon1, epsilon2)
for(t in 3:Tmax) {
    x[t] <- a %*% x[t-1:2] + epsilon[t] # %*% はベクトルの内積計算
}
x |> ts() |> as_tsibble() |> autoplot(value)

複数の系列を表示するには，以下に一連の手続きを記述して関数化しておくのが良い．

my_ar <- function(a, epsilon){ 
    p <- length(a) # 次数pを取得
    Tmax <- length(epsilon) # 時系列の長さを取得
    x <- double(Tmax)      # 変数を用意
    x[1:p] <- epsilon[1:p] # 初期値は(epsilon1,...)
    for(t in (p+1):Tmax) {
        x[t] <- a %*% x[t-1:p] + epsilon[t]
    }
    return(x) # 計算結果のxを返す
}
#' 使い方は a と epsilon(ホワイトノイズ)を指定する
x <- my_ar(a = c(0.6, 0.3, 0.1), epsilon = rnorm(100))
x |> ts() |> as_tsibble() |> autoplot(value)

関数の引数として Tmax を指定する方法もあるが，様々な分布のホワイトノイズを試したい場合もあるので生成の元となるノイズを渡す形で定義してある．

データフレームを作成して表示する．後の演習で作成した時系列データを利用するので，適当なオブジェクトに保存しておく．

ts_ar <- ts(replicate(K, my_ar(a = a, epsilon = rnorm(Tmax))))
ts_ar |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    labs(title = "AR(2)", x = "Time", y = "Observations")

MA(2)モデルのシミュレーションを行う．

b <- c(0.44, 0.08) # MAの係数
epsilon <- rnorm(Tmax)
x <- ts(double(Tmax))
x[1:2] <- epsilon[1:2]   
for(t in 3:Tmax) {
    x[t] <- b %*% epsilon[t-1:2] + epsilon[t]
}
x |> ts() |> as_tsibble() |> autoplot(value)

同様に複数の系列を表示する．

my_ma <- function(b, epsilon){
    q <- length(b) # 次数qを取得
    Tmax <- length(epsilon) # 時系列の長さを取得
    x <- double(Tmax)
    x[1:q] <- epsilon[1:q]   
    for(t in (q+1):Tmax) {
        x[t] <- b %*% epsilon[t-1:q] + epsilon[t]
    }
    return(x)
}
ts_ma <- ts(replicate(K, my_ma(b = b, epsilon = rnorm(Tmax))))
ts_ma |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    labs(title = "MA(2)", x = "Time", y = "Observations")

ARMA(2,1)モデルのシミュレーションを行う．

a <- c(0.8, -0.64) # ARの係数
b <- -0.5          # MAの係数
epsilon <- rnorm(Tmax)
x <- double(Tmax)
x[1:2] <- epsilon[1:2]
for(t in 3:Tmax) {
    x[t] <- a %*% x[t-1:2] + b %*% epsilon[t-1] + epsilon[t]
    #' bは1次元なのでこの問題では b*epsilon でも可
}
x |> ts() |> as_tsibble() |> autoplot(value)

複数の系列を表示する．

my_arma <- function(a, b, epsilon){
    p <- length(a)
    q <- length(b)
    r <- max(p,q)
    Tmax <- length(epsilon) # 時系列の長さを取得
    x <- double(Tmax)
    x[1:r] <- epsilon[1:r]
    for(t in (r+1):Tmax) {
        x[t] <- a %*% x[t-1:p] + b %*% epsilon[t-1:q] + epsilon[t]
    }
    return(x)
}
ts_arma <- ts(replicate(K, my_arma(a = a, b = b, epsilon = rnorm(Tmax))))
ts_arma |> as_tsibble() |>
    autoplot(value, show.legend = FALSE) +
    labs(title = "ARMA(2,1)", x = "Time", y = "Observations")

ここでは時系列の生成過程を知ってもらうために自作の関数を作成したが，関数 stats::arima.sim() や stats::filter() などを利用することもできる．上記のARMA(2,1)のシミュレーションは以下のようにして行うことができる．

arima.sim(model = list(ar = c(0.8, -0.64), # ARの係数ベクトル
                       ma = c(-0.5)),      # MAの係数ベクトル
          n = Tmax, # 時系列の長さ
          innov = epsilon) # 乱数系列を渡す場合(渡さなければ標準正規乱数が使われる)

詳細は ?stats::arima.sim() を参照のこと

自己相関

自己相関や相互相関および共分散の計算には関数 feats::ACF() を用いればよい．

相関

ACF(.data,  y, ...,  lag_max = NULL,
    type = c("correlation", "covariance", "partial"),
    na.action = na.contiguous, demean = TRUE, tapered = FALSE)
#' .data: 時系列データ (tsibbleクラス)
#' y: 計算対象の列名
#' type: 標準は相関, 共分散と偏相関を選ぶこともできる
#' na.action: 欠損値の処理，標準は欠損を含むと計算しない
#' demean: 共分散の計算において平均を引くか否か
#' 詳細は '?feats::ACF' を参照

関数 ACF()の返値を関数 autoplot() に渡せば適切なグラフを描画することができる．

toy_acf <- arima.sim(model = list(ar = c(0.8, -0.64),
                                  ma = c(-0.5)),
                     n = 200) |>
    as_tsibble() |> ACF(value) 
toy_acf |> autoplot()

base R にも自己相関を計算する関数は用意されている．

acf(x, lag.max = NULL,
    type = c("correlation", "covariance", "partial"),
    plot = TRUE, na.action = na.fail, demean = TRUE, ...)
#' x: 時系列データ
#' lag.max: 計算するラグの最大値
#' type: 標準は相関, 共分散と偏相関を選ぶこともできる
#' plot: 描画するか否か
#' na.action: 欠損値の処理，標準は欠損を含むと計算しない
#' demean: 共分散の計算において平均を引くか否か
#' 詳細は '?stats::acf' を参照

引数 plot の真偽で描画するか，計算のみするかを制御できる．

問題

以下の問に答えなさい．

同じAR過程のモデルから生成した時系列の自己相関を比較しなさい．(前の練習問題を利用すればよい)
MA過程についても同様な比較を行いなさい．
ARMA過程についても同様な比較を行いなさい．

解答欄

解答例

格子状に表示する時系列の数を設定する．

K <- 4 # 4つを並べて比較する

前問で作成したAR(2)モデルの自己相関を計算し，視覚化する．

ts_ar |> as_tsibble() |> ACF(value) # 計算結果はtsibbleクラスになる

# A tsibble: 115 x 3 [1]
# Key:       key [5]
   key           lag   acf
   <chr>    <cf_lag> <dbl>
 1 Series 1        1 0.844
 2 Series 1        2 0.783
 3 Series 1        3 0.654
 4 Series 1        4 0.564
 5 Series 1        5 0.466
 6 Series 1        6 0.388
 7 Series 1        7 0.314
 8 Series 1        8 0.256
 9 Series 1        9 0.193
10 Series 1       10 0.151
# ℹ 105 more rows

ts_ar |> as_tsibble() |> ACF(value) |> autoplot()

ts_ar[,1:K] |> as_tsibble() |> ACF(value) |>
    autoplot() + # 格子状に並べるには facet を指定する
    facet_wrap(key ~ ., nrow = 2) +
    labs(title = "AR(2)")

tsibble クラスのデータフレームを保存しておいて利用することもできる．その場合は関数 filter などを使って操作することができる．

tsb_ar <- ts_ar |> as_tsibble()
tsb_ar |> ACF(value) |> autoplot()

tsb_ar |>
    filter(str_detect(key,"[1-4]")) |> # Series [1-4] のみ取り出す
    ACF(value) |>
    autoplot() + 
    facet_wrap(key ~ ., nrow = 2) +
    labs(title = "AR(2)")

同様にMA(2)モデルの自己相関を計算し，視覚化する．

ts_ma[,1:K] |> as_tsibble() |> ACF(value) |>
    autoplot() + 
    facet_wrap(key ~ ., nrow = 2) +
    labs(title = "MA(2)")

ARMA(2,1)モデルの自己相関を計算し，視覚化する．

ts_arma[,1:K] |> as_tsibble() |> ACF(value) |>
    autoplot() + 
    facet_wrap(key ~ ., nrow = 2) +
    labs(title = "ARMA(2,1)")