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

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

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

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

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

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

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

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

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

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

(2)　時系列の図示

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

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

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

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

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

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

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

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

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

(4)　差分

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

図4　UKgas の差分のプロット

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

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

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

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

で定義され、標本自己相関関数は

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

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

図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 のトレンドが直線ではなく、指数分布に近似してるからである。

　

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

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

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

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

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

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

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

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

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