[連載]フリーソフトによるデータ解析・マイニング第33回

Rと時系列(1)

1.時系列とは

 時間とともに変動する現象に対して時間の順序で測定・観測した結果の記録を時系列データと言い、略して時系列(time series)と言う。時系列データは多くの分野で様々な目的で取り扱われる。日常の社会生活の中でよく見受けられるものには、心電図や脳波のような医療データ、気温や気圧のような気象データ、株価および為替レートのような金融・経済データなどがある。
 時系列データは、常に変動を伴うものである。その振る舞いを統計的に分析し、データ変動の特徴を捉え、現象の解明と将来の変動を予測・制御しようとするのが時系列データ分析の主要な目的である。
 ちなみに、2003年ノーベル経済学賞の受賞の対象となった内容は、経済時系列分析に関するものである。

2.時系列データの形式と操作

 時系列データは,通常ある一定の離散的時間間隔で測定・観測したものである。観測・測定値 に関する標本サイズ の時系列データは、測定した時間 をインデックスに次のように表記される。

y1,y2,…,yt-k,…,yt-1,yt,yt+1,…,yt+k,…,yn-1,yn

(1)Rにおける時系列データの属性と表示

 Rでは、マトリックス、データフレームに並んで時系列データオブジェクト形式がある。
 基本的な時系列データの操作と分析に関する関数はパッケージstatsに用意されている。パッケージstatsはRを起動して直接に利用することが可能となっている。
 本稿では、Rの中に用意されている時系列データを用いて説明することにする。Rには、ある女性の血液中の黄体形成ホルモンについて10分の間隔で測定した48のデータlhがある。次のコマンドが返す結果からデータlhの属性が確認できる。返された結果tsは時系列(time series)の頭文字で、属性が時系列であることを意味する。

> class(lh)
[1] "ts"

 Rにおける時系列データオブジェクトtsは、データに時間などの情報を加えて構造化したものである。例えば、次に示すデータlhの内部には、時系列が始まる時間を1、終わる時間を48として、観測データが記録されている。

>lh
lhの結果

 データlhでは10分を1つの時間間隔としているが、この時間内では、1回しか観測していないのでfrequency =1になっている。1年を単位とし、月ごとに観測したデータの場合はfrequency =12、四半期ごとに観測した場合は、frequency =4となる。
 Rには1960年から1986年までのイギリスのガス消費量を四半期ごとに観測した時系列データUKgasがある。関数start、 endを用いて時系列データの開始時間と終了時間を、関数frequencyで時間ごとに観測した回数を返すことができる。

>start(UKgas)
[1]1960   1
>end(UKgas)
[1]1986   4
>frequency(UKgas)
[1]4
>UKgas

  Qtr1 Qtr2 Qtr3 Qtr4
1960 160.1 129.7 84.8 120.1
1961 160.1 124.9 84.8 116.9
<中略>        
1985 1087.0 534.7 281.8 787.6
1986 1163.9 613.1 347.4 782.8

 Rでは関数windowを用いて、時系列の一部分を切り出すことができる。次にデータUKgasを用いた関数windowの使用例を示す。

> window(UKgas,c(1975,2),c(1979,3))

  Qtr1 Qtr2 Qtr3 Qtr4
1975   321.8 177.7 409.8
1976 593.9 329.8 176.1 483.5
1977 584.3 395.4 187.3 485.1
1978 669.2 421.0 216.1 509.1
1979 827.7 467.5 209.7  

(2)時系列の図示

 時系列データは関数ts.plotを用いて横軸を時間、縦軸を観測・測定値とした折れ線グラフを作成することができる。関数ts.plotは、plotに略して用いることもできる。関数plotは自動的にオブジェクトの属性を識別することができるので、オブジェクトが時系列の場合は自動的に関数ts.plotの機能を用いて作図される。次にデータlh、UKgasの時系列グラフの作成例を示し、その結果を図1、2に示す。

>ts.plot(lh)


図1 データlhの時系列プロット

図1 データlhの時系列プロット

>ts.plot(UKgas)


図2 データUKgasの時系列プロット

図2 データUKgasの時系列プロット

 関数ts.plotでは、通常の作図関数plotと同じく線種、色などを自由に設定することができる。
 Rには、1974年から1979年までのイギリスで喘息、気管支炎、肺気腫による死亡数を月ごとに記録した時系列データldeaths、これを男女別に分けたmdeaths、fdeathsがある。この3種類のデータを1つのグラフに表すコマンドを次に示し、その結果を図3に示す。

> ts.plot(ldeaths, mdeaths, fdeaths,gpars=list(xlab="年", ylab="死亡数",lty=c(1:3),col=c(1:3)))
> legend(locator(1),c("全体","男性","女性"),lty=c(1:3),col=c(1:3))

 関数legendの中の引数locator(1)を用いることにより、作成した画面に凡例を置く場所をマウスで指定することができる。locator(2)にすると凡例の外枠が描かれない。


図3 3種類のdeathsの時系列プロット

図3 3種類のdeathsの時系列プロット

(3)データオブジェクトの作成

 非時系列属性のデータは、関数ts、as.tsを用いて時系列属性に変換することができる。関数tsを用いる際は、開始時間(start)、時間単位ごとの観測数(frequency)を引数で指定することが必要である。例えば、1から120までの整数を1995年から2004年まで、1年に12回観測したデータとして時系列オブジェクトを作成するには、次のように関数tsを用いる。

>temp<-ts(1:120,start=c(1995,6),frequency=12)
> class(temp)
[1]"ts"
>temp
tempの結果

 上記のコマンドの例でわかるように、引数startは時系列が始まる時間を指定する。用いる全てのデータを時系列に変換するときには、引数startを指定すると終わる時間は唯一であるので、引数endを指定しなくてもよい。
 関数as.tsはデータマトリックスやデータフレームを時系列の属性に変換するのに便利である。

(4)ラグ処理

 時系列データを分析する際、データの時間をシフトして比較したりすることも少なくない。時間の遅れのことを時系列ではラグ(lag)と呼ぶ。
 時間tを基準としたとき、yt-1ytの1次ラグ、yt-2を2次ラグ、yt-kk次ラグと呼ぶ。
 Rでは、ラグ処理を行う関数lag(x,k=)がある。引数kは正負の整数である。次に、データldeathsを用いたラグ処理の結果を示す。

>ldeaths
ldeathsの結果
>lag(ldeaths,k=5)
lag(ldeaths,k=5)の結果

 また、Rには時系列データの併合などを行う関数ts.union、ts.intersectが用意されている。

(4)差分

 値ytからyt-1を引いた値⊿ytytyt-1を差分(階差)と呼ぶ。時系列データにトレンドを含む場合は、差分操作で線形関係のトレンドを除去することができる。
 Rには差分を求める関数diffが用意されている。次に時系列データUKgasの差分を図示するコマンドを示し、その結果を図4に示す。図2と比べると分かるように右肩上がりの時系列データの差分⊿ytではそのトレンドが除去されている。

> plot(diff(UKgas))


図4 UKgasの差分のプロット

図4 UKgasの差分のプロット

 関数diffには、引数lagがあり、デフォルトは1になっている。このlagは、ラグ関数のlagと同じ意味を持っている。関数diffの引数lagと関数lagの引数kとの対応関係についてデータUKgasを用いて次に示す。

diff(UKgas,lag=2)=UKgas-lag(UKgas,k=-2)

 関数diffには、差分演算を繰り返す引数differencesがある。デフォルトは1になっている。

3.自己共分散と自己相関

 通常のデータの基本的な統計的特性を表す統計量として、平均、分散、相関などが用いられている。これに対応して時系列では、平均、自己共分散(autocovariance)、自己相関(autocorrelation)と呼ばれる統計量がある。
 時系列{y1,y2,…,yn-1,yn}の標本平均は

時系列{y1,y2,…,yn-1,yn}の標本平均


で定義される。時系列の平均や自己共分散の性質が、時間が経過しても変化しない場合、その時系列を定常時系列と言い、それらの性質が何らかの形で変化するものを非定常時系列と言う。
 定常時系列について標本データのytyt-kとの標本自己共分散関数は

標本自己共分散関数


標本自己相関関数は

標本自己相関関数


で定義される。
 Rには自己共分散、自己相関を求める関数acfがある。

acf(x, type =””, plot = TRUE,...)

 引数xは時系列データ、引数typeには"correlation"(自己相関), "covariance"(自己共分散), "partial"(偏相関)の中から1つを選んで指定することができる。デフォルトは" correlation "になっている。引数plotは図示の指定引数で、デフォルトはplot = TRUEになっているので、自動的に作図される。図示のデフォルトは関数plot.acfのデフォルト値に従う。また、欠損値の扱いに関する引数na.actionなどがある。
 データUKgasを用いた関数acfの使用例を次に示し、その結果を図5に示す。

>acf(UKgas)

図5 UKgasの自己相関プロット7

図5 UKgasの自己相関プロット7

 図5から、横軸1、2、3、4、5のところに自己相関係数の値が大きく、データUKgasは1年単位の周期性が顕著であることが分かる。
 図5の横の破線は95%の信頼区間である。信頼区間は引数ciで指定する。例えば、信頼区間を90%にするときはacf(UKgas,ci=0.9)のように書く。
 データにおける周期性とトレンドは関数diffを用いて除去することも可能である。 関数diffの引数lagを1、2、3、4にしたデータUKgasの差分を求めるコマンドを次に示し、その結果を図6に示す。ここではまずUKgasを対数変換し、差分を計算している。これは図2で分かるようにデータUKgasのトレンドが直線ではなく、指数分布に近似してるからである。

> par(mai=rep(0.2,4),mfrow=c(4,1))
> for(i in 1:4)plot(diff(log(UKgas),lag=i))

図6 異なる差分のプロット

図6 異なる差分のプロット

 データUKgasは1年を四半期ごとに記録したデータであり、年周期を持っているので、関数diffの引数lagを4にした差分には周期成分が除去される。
 また異なる2つの時系列データ間の相互共分散および相互相関は関数ccfで求めることができる。次にデータmdeathsとfdeathsとの相互相関を求める関数ccfの使用例を示し、その結果を図7に示す。

> ccf(mdeaths, fdeaths)

図7 mdeathsとfdeathsの相互相関プロット

図7 mdeathsとfdeathsの相互相関プロット

 図7から両時系列データは強い相関と周期性を持っていることが読み取れる。
 偏自己共分散・相関は関数pacfで求めることができる。

4.スペクトル分析

 時系列データは、幾つかの成分が混合されているものと考えられる。例えば、トレンド、周期の成分、ノイズなど。時系列の分析では、周期成分の有無やその成分の値などを分析する必要がある。時系列データに隠されている周期性を解析する方法としてスペクトル解析がある。
 時系列における自己共分散Ckのフーリエ変換(Fourier transform)が可能であるとき、周波数-1/2≦f≦1/2上で定義される関数

関数

をパワースペクトル密度関数(power spectral density function)、略してスペクトル(spectrum)と呼ぶ。また、次式のように上記のパワースペクトル密度関数におけるCkに代わって、標本データy1,y2,…,ynの標本自己共分散標本自己共分散を用いて定義したものをピリオドグラム(preiodogram)と呼ぶ。

関数

 ただし、周波数は 、 である。ピリオドグラムはスペクトル密度を推定する1つのツールとして用いられている。
 Rには、ピリオドグラムを用いてスペクトル密度関数を推定する関数spec.pgramがある。関数spec.pgramは、FFT(fast Fourier transform)でピリオドグラム(periodogram)を求め、対数軸(縦軸)のグラフを作成し、引数spansを用いて修正Daniell平滑化を行うことが可能でる。
 次に関数spec.pgramによるデータUKgas、ldeathsのピリオドグラムと1組のパラメータによる平滑化のコマンドを示し、その結果を図8に示す。

par(mfrow=c(2,2))
spec.pgram(UKgas)
spec.pgram(UKgas,spans=c(3,3))
spec.pgram(ldeaths)
spec.pgram(ldeaths,spans=c(3,3))


図8 ピリオドグラムのプロット

図8 SVMのイメージ

 関数spec.pgramによるピリオドグラムのプロットでは、右側の縦のバーは95%信頼区間を示すものである。
 引数spansに用いるパラメータは移動平均の長さで、奇数にする必要がある。スペクトルに関するより詳しい参考文献としては[1]と[2]がある。
 Rにはスペクトルを求める関数spectrumがある。この関数spectrumは、関数spec.pgramとspec.arを統合したもので、デフォルトにはspec.pgramが設定されている(method="pgram")。関数spec.arは自己回帰によるスペクトル解析の関数である。次にデータUKgas、ldeathsを用いた自己回帰による関数spectrumの使用例を示し、その結果のプロットを図9に示す。

> par(mfrow=c(1,2))
> spectrum(UKgas,method="ar")
> spectrum(ldeaths,method="ar")

図9 自己回帰のスペクトルプロット

図9 自己回帰のスペクトルプロット



参考文献
[1] 伊藤幹夫、大津泰介、戸瀬信之、中東雅樹 訳:S-PLUS による統計解析、シュプリンガー・フェアラーク東京
[2] 北川源四郎:時系列解析入門、岩波書店