Chap4 時系列分析

コンテンツ

4.1 時系列分析の概要

S-PLUSの組み込みのデータセットである．"co2"を使って，時系列分析とは何かを 具体的に説明する．

"co2"は米国のオークリッジ国立研究所が保持している気象データベースからとったもので，1959年1月から1990年12月までの二酸化炭素の濃度のデータである．

"co2"の内容と，概形は下の通りである．

         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct    Nov    Dec
1959: 315.42 316.32 316.49 317.56 318.13 318.00 316.39 314.66 313.68 313.18 314.66 315.43
1960: 316.27 316.81 317.42 318.87 319.87 319.43 318.01 315.75 314.00 313.68 314.84 316.03
1961: 316.73 317.54 318.38 319.31 320.42 319.61 318.42 316.64 314.83 315.15 315.95 316.85
1962: 317.78 318.40 319.53 320.41 320.85 320.45 319.44 317.25 316.12 315.27 316.53 317.53
1963: 318.58 318.92 319.70 321.22 322.08 321.31 319.58 317.61 316.05 315.83 316.91 318.20
1964: 319.41 320.07 320.74 321.40 322.06 321.73 320.27 318.54 316.54 316.71 317.53 318.55
1965: 319.27 320.28 320.73 321.97 322.00 321.71 321.05 318.71 317.65 317.14 318.71 319.25
1966: 320.46 321.43 322.22 323.54 323.91 323.59 322.26 320.21 318.48 317.94 319.63 320.87
1967: 322.17 322.34 322.88 324.25 324.83 323.93 322.39 320.76 319.10 319.23 320.56 321.80
1968: 322.40 322.99 323.73 324.86 325.41 325.19 323.97 321.92 320.10 319.96 320.97 322.48
1969: 323.52 323.89 325.04 326.01 326.67 325.96 325.13 322.90 321.61 321.01 322.08 323.37
1970: 324.34 325.30 326.29 327.54 327.54 327.21 325.98 324.42 322.91 322.90 323.85 324.96
1971: 326.01 326.51 327.01 327.62 328.76 328.40 327.20 325.28 323.20 323.40 324.64 325.85
1972: 326.60 327.47 327.58 329.56 329.90 328.92 327.89 326.17 324.68 325.04 326.34 327.39
1973: 328.37 329.40 330.14 331.33 332.31 331.90 330.70 329.15 327.34 327.02 327.99 328.48
1974: 329.18 330.55 331.32 332.48 332.92 332.08 331.02 329.24 327.28 327.21 328.29 329.41
1975: 330.23 331.24 331.87 333.14 333.80 333.42 331.73 329.90 328.40 328.17 329.32 330.59
1976: 331.58 332.39 333.33 334.41 334.71 334.17 332.88 330.77 329.14 328.77 330.14 331.52
1977: 332.75 333.25 334.53 335.90 336.57 336.10 334.76 332.59 331.41 330.98 332.24 333.68
1978: 334.80 335.22 336.47 337.59 337.84 337.72 336.37 334.51 332.60 332.37 333.75 334.79
1979: 336.05 336.59 337.79 338.71 339.30 339.12 337.56 335.92 333.74 333.70 335.13 336.56
1980: 337.84 338.19 339.90 340.60 341.29 341.00 339.39 337.43 335.72 335.84 336.93 338.04
1981: 339.06 340.30 341.21 342.33 342.74 342.07 340.32 338.27 336.52 336.68 338.19 339.44
1982: 340.57 341.44 342.53 343.39 343.96 343.18 341.88 339.65 337.80 337.69 339.09 340.32
1983: 341.20 342.35 342.93 344.77 345.58 345.14 343.81 342.22 339.69 339.82 340.98 342.82
1984: 343.52 344.33 345.11 346.88 347.25 346.61 345.22 343.11 340.90 341.17 342.80 344.04
1985: 344.79 345.82 347.25 348.17 348.75 348.07 346.38 344.52 342.92 342.63 344.06 345.38
1986: 346.12 346.79 347.69 349.38 350.04 349.38 347.78 345.75 344.70 344.01 345.50 346.75
1987: 347.86 348.32 349.26 350.84 351.70 351.11 349.37 347.97 346.31 346.22 347.68 348.82
1988: 350.29 351.58 352.08 353.45 354.08 353.66 352.25 350.30 348.58 348.74 349.93 351.21
1989: 352.62 352.93 353.54 355.27 355.52 354.97 353.74 351.51 349.63 349.82 351.12 352.35
1990: 353.47 354.51 355.18 355.98 356.94 355.99 354.58 352.68 350.72 350.92 352.55 353.91

図x.co2の概形

時系列分析の第一歩は，概形を眺めてデータの挙動を把握することである． 図xは下の2つの変動に分けることができる．

1. 二酸化炭素の濃度は全体的にみると増加している
2. 二酸化炭素の濃度は1年周期で細かい増減を繰り返している

前者をトレンド成分，後者を季節成分と呼ぶ． 時系列分析はこの2つの変動の影響を取り除き， ノイズ成分を抽出してモデルを構築して， 最後にトレンド成分および季節成分を再び考慮することで行われる．

トレンド成分および季節成分を取り除いたノイズ成分が， 定常であれば，後で説明する時系列モデルを適用することができ， 構築したモデルを用いて， 観測期間より後の二酸化炭素濃度の予測をすることが可能となる．

4.2 時系列分析のための基礎知識

ここでは，定常という概念と，分析に必要な基礎知識について説明する．

時系列がある． 次の2つの条件を満たすときに，は定常であるという．

まず平均を次のように定義する．

このときはすべてのに対して一定となる．

またとの共分散を次のように定義する．

このときは各について，すべてのに対して一定となる．

つまり時系列において，いかなるときでも平均と共分散が等しければ，は定常であるといえる．

を定常時系列とする． ここで，今後の分析に必要となる2つの関数を定義する．

ラグにおける， 自己共分散関数(ACVF:AutoCoVariance Function) は次のように定義する．

ラグにおける，自己相関係数(ACF:AutoCorrelation Function)は次のように定義する．

4.3 時系列モデル

時系列分析をする大きな目的のひとつが，将来の予測をすることである． 時系列モデルによって，与えられた時系列の特徴を捕らえることができれば， そのモデルによって，多くの将来のシナリオを作ることもできるし， ある信頼区間内で，将来どのような挙動をするか予測することもできる． 以下では，代表的な時系列モデルを紹介する．

4.3.1 IIDノイズ

一番基本的な時系列モデルとして， IID(Independently and Identically Distributed)ノイズがある． IIDノイズとは，観測値が独立同一分布に従う，平均のモデルである． このモデルは直接予測はできないが，他のモデルを構築する上で基礎となる．

IIDノイズのACVFは次の通りである．

平均が，ACVFが時間によって変化しないことから， IIDノイズは定常であるといえる．

が平均，分散がのIIDノイズであることを次のように表現する．

4.3.2 ホワイトノイズ

ホワイトノイズは観測値が無相関な分布に従う，平均のモデルである．

ホワイトノイズのACVFはIIDノイズと同じで， 平均が，ACVFが時間によって変化しないことから， ホワイトノイズも定常であるといえる．

が平均， 分散がのホワイトノイズであることを次のように表現する．

ここでホワイトノイズの定常性について確認する． S-PLUSでは"rnorm"という関数によって，平均， 分散に従う， ホワイトノイズを発生することができる． 500のホワイトノイズを発生して，"wn"という変数に付値する．

wn = rnorm(500)

"wn"を用意したら，自己相関をプロットするメニューを呼び出す． メニューバーから，[統計]→[時系列解析]→[自己共分散・相関]の順番で選択する．

すると，以下のようなダイアログが現れる． ここで，どのデータの自己共分散あるいは自己相関を見るのかを設定する．

このウィンドウで，次のような設定を行う．

• "wn"というデータを使用する(Data,Object / Variable)
• 自己共分散か自己相関の選択(Options,Estimate Type)
• 結果の出力(Results,Plot Results)

はじめに，[Object]には使用するデータセットを入力（選択）する． 例では"wn"と入力する． [Variable]は使用するデータセットが複数の列からなる場合， 必要な変数列を入力（選択）する． 例では"wn"は1変数のみのデータセットなので， 自動的に"wn"が選択される．

次に[Estimate Type]で自己共分散か自己相関か出力したいものを選択する． 例では自己相関が見たいので，"autocorrelation"を選択する．

最後に[Plot Results]のチェックボックスをオンにすると， 先に選択した自己共分散あるいは自己相関をグラフ化する． 例では，[Plot Results]のチェックボックスをオンにしている．

図x.ホワイトノイズの自己相関

図xの横軸はlag，縦軸は自己相関の大きさである．また2本の点線は95%信頼区間である． lag0の自己相関は必ず1になるが，その他のlagにおいて，自己相関の絶対値はすべて信頼区間内に収まっている．これは，"wn"という系列が，定常であることを支持している． なお，自己相関の絶対値が2，3個程度信頼区間外に落ちた場合も定常であると考えられる．

4.3.3 差分モデル

時系列を時系列で引くことを時系列の（1階）差分をとるという． を回，差分を繰り返してとったものを，の階差分と呼びと表す． 階の差分をとると，個のデータが失われる．

S-PLUSでは差分をとる関数"diff"が用意されている． この関数を用いて実際に"co2"の差分をとってみる． まず"co2"の1階差分をとったものを"co2.diff"という変数に付値する．

> co2.diff = diff(co2)

"length"という関数を用いて，"co2"と"co2.diff"のデータ数を確認する．

> length(co2)
[1] 384
> length(co2.diff)
[1] 383

次に"co2.diff"の概形を確認する．

図x.co2の1階差分

"co2"の増加傾向のトレンドが取り除かれていることがわかる． このように差分をとることで，データのトレンドを除去することができる．

4.3.4 AR(AutoRegressive:自己回帰)モデル

が， 次の式で表されるとき， 過程であるという．

4.3.5 MA(Moving Average:移動平均)モデル

が， 次の式で表されるとき， 過程であるという．

4.3.6 ARMA(AutoRegressive Moving Average:自己回帰移動平均)モデル

が， 次の式で表されるとき， 過程であるという．

4.3.7 ARIMA(AR Integrated MA:自己回帰和分移動平均)モデル

階の差分をとることで， がとなるような過程を， 過程であるという．

4.4 ARIMAモデルの構築

今回用いるのは，架空のアイスクリーム店のデータである．この店は1991年の1月にオープンし，オープン時から毎月の来客数を記録している．その一部が以下の表xであるものとする．（データが欲しい方は「こちら(icecream.csv)」からどうぞ）

データをS-PLUSにインポートしなくてはなりませんが，データのインポート・エクスポートをご参照ください．なお，その際データセット名を"icecream"にしていただきますと，今後データセット名を気にすることなく進めていくことができます．

S-PLUSにおいて時系列解析を行うためには，データセットを時系列のクラスにしておく必要がある．1991年1月から始まる月次の時系列データに変換するには，コマンドウィンドウに下のように入力して，実行する．

> icecream = rts(icecream, start=1991, frequency=12, units="months")

rtsは一定間隔の時系列を生成する関数であり，startで1991を指定することにより，1991年から始まる時系列データとなる．また，frequencyとunitsにそれぞれ，12と"months"を指定することにより，月次の時系列データとなる．

次にコマンドウィンドウに下のように入力し実行する．

> icecream

すると下のような出力が返ってきて，時系列データが作成できたことが確認できる．

       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1991: 3296 2882 2860 3071 3351 3980 4014 4008 3651 3415 3320 3092
1992: 2979 3140 3156 3004 4334 3967 4264 3807 4627 3616 3087 2892
1993: 3434 3235 3806 3206 4461 4259 4502 4356 4350 3917 3688 3180
1994: 3636 3936 3687 3720 4196 4467 4539 4735 4326 4238 3779 3440
1995: 3924 3489 3915 4269 4624 4752 4792 4718 4801 4430 4092 4158
1996: 4058 3951 4086 4317 5213 5148 5125 5197 4963 4652 4092 4185
1997: 4266 4336 4243 4203 5195 4719 5098 5561 5292 4785 4417 4554
1998: 4645 4623 4212 4220 5043 5545 5544 5630 4912 4996 4522 4526
1999: 4520 4390 4645 4828 5503 5273 5881 5356 5475 5311 4801 4925
2000: 4871 4869 4720 4991 5537 5622 5828 5627 5742 5804 4966 5297
2001: 4684 4997 5163 5246 5553 6025 6017 5841 5954 5756 5642 4990
2002: 5328 4930 5262 5545 5762 6041 6613 6516 5987 5934 5531 5704
2003: 5163 5225 5171 5468 5895 6677 6484 6550 6303 6076 5542 5483
2004: 5742 5543 5904 5867 6084 6712 7062 6496 6324 6060 6165 5778
2005: 5980 5624 5962 6072 6342 6623 6946 6807 6690 6247 6146 6432
 start     deltat frequncy 
  1991 0.08333333        12
 Time units :  months

さらにコマンドウィンドウに下のように入力し実行すると，このデータのグラフを描くことができる．"tsplot"は時系列データのグラフを作図する関数である．図xは出力された"icecream"の概形である．

> tsplot(icecream)

図x.icecreamの概形

時系列分析はデータを視覚化して，把握することが重要である． 図xは下の2つの変動に分けることができる．

1. 来客数は全体的にみると増加している
2. 来客数は夏が多く，冬が少ないというように1年周期で細かい増減を繰り返している

前者をトレンド成分，後者を季節成分と呼ぶ． ARIMAモデルを用いた時系列分析はこの2つの変動の影響を取り除き， 定常なノイズ成分を抽出してモデルを構築して， 最後にトレンド成分および季節成分を再び考慮することで行われる．

データの概形を確認したら，ARIMAモデルを用いた時系列分析を行うメニューを呼び出す． メニューバーから，[統計]→[時系列解析]→[ARIMAモデル]の順番で選択する．

図x.ARIMAモデルを選択

すると，以下のようなダイアログが現れる． まず[Model]タブをクリックする．ここで，どのようなARIMAモデルを構築するのか設定する．

図x.ARIMAモデルウィンドウ(Model)

このウィンドウはARIMAモデルを構築する際の設定メニューである． これから，次のような設定を行う．

• "icecream"というデータを使用する(Data,Object / Variable)
• ARIMAモデルの次数，階数(ARIMA Model Order,Autoregressive / Difference / Moving Ave.)
• 季節成分を考慮するか(ARIMA Model Periodicity,Seasonality / Period)
• 結果の出力をする(Results,Save As / Print Results)

はじめに，[Object]には使用するデータセットを入力（選択）する． 例では"icecream"と入力する．[Variable]は使用するデータセットが複数の列からなる場合， 必要な変数列を入力（選択）する．例では"icecream"は1変数のみのデータセットなので， 自動的に"icecream"が選択される．

次にデータセットをあてはめるARIMAモデルの次数と階数を入力（選択）する． [Autoregressive]は自己回帰の次数，[Difference]は差分をとる階数， [Moving Avg.]は移動平均の次数である．例では[Autoregressive]を1， [Difference]を1，[Moving Avg.]を1と入力している． これは，ARIMA(1,1,1)モデルを構築することを意味している．

さらに季節成分を考慮するかを選択する．[None(1)]は季節成分を考慮しない， [Quarterly(4)]は季節成分の周期を四半期（3ヶ月）として考慮する， [Monthly(12)]は季節成分の周期を1年として考慮する， そしてOtherは季節成分を任意の周期で考慮する． 例では季節成分を考慮しない[None(1)]を選択している． [Period]は[Seasonality]で[Other]を選択した場合， どのくらいの期間を周期とするかを入力（選択）できる．

最後に[Save As]は出力結果を保存したいときに入力する． 入力した名前で，arimaというクラスのobjectを得ることができる． 例では，実行すれば"icecream.results"というobjectを得る． [Print Results]のチェックボックスをオンにすると出力結果をReportとして表示する． 例ではチェックボックスをオンにしている．

モデルの構築に関して設定が終わったら， 次に[Diagonostics]タブをクリックする． ここでモデルの評価のためにグラフの出力を設定する．

図x.ARIMAモデルウィンドウ(Diagonostics)

このウィンドウはモデルの評価のためのグラフの出力をする際の設定メニューである． これから，次のような設定を行う．

• 残差の自己相関関数（以下ACF:AutoCorrelation Function）と偏自己相関関数（以下PACF:Partial ACF）の値を出力をする(Save,Save As/Autocorrelation of Residuals)
• 残差のACFとPACFをグラフとして出力する(Plots,Plot Diagnostics)

まず[Save As]は残差のACFとPACFのグラフを 再現するための値を保存したいときに入力する． 入力した名前で，namedというクラスのobjectを得ることができる． 例では，実行すれば"icecream.diagnostics"というobjectを得る． また，ACFとPACFの値を出力するために， [Autocorrelation of Residuals]のチェックボックスをオンにする． 例では，チェックボックスをオンにしている．

次に[Plot Diagnostics]のチェックボックスをオンにすると 先に出力したACFとPACFの値をグラフ化する． 例では，[Plot Diagnostics]のチェックボックスをオンにしている．

モデルの評価に関して設定が終わったら， 次に[Forecast]タブをクリックする． ここで将来予測のグラフの出力を設定する．

図x.ARIMAモデルウィンドウ(Forecast)

このウィンドウは将来予測のグラフの出力をする際の設定メニューである． これから，次のような設定を行う．

• 予測期間(Time Periods To Forecast)
• 予測期間の推移を出力する(Plots)
• 予測結果を保存する(Save Forecast As)

まず[Time Periods To Forecast]には，予測期間を入力する． 例では，24としているので24ヶ月（2年）の予測をする．

次に，[Plot Forecast]のチェックボックスをオンにすると出力結果をグラフで出力する．． 例ではチェックボックスをオンにしている．

最後に[Save Forecast As]は予測を再現するための値を保存したいときに入力する． 入力した名前で，data.frameというクラスのobjectを得ることができる． 例では，実行すれば"icecream.forecast"というobjectを得る．

これで設定はすべて終了である．画面左下の[OK]をクリックする．

4.4 ARIMAモデルを用いた時系列分析

設定終了後，[OK]をクリックすると，2つのウィンドウが新しく立ち上がる．

1. ARIMAモデルの統計量に関するレポート
2. 残差のACFとPACFのグラフおよび予測のグラフ

4.4.1 レポートの読み方

まず，ARIMAモデルの統計量に関するレポートについて核となる部分を解説する． 以下のように記述されているウィンドウをクリックする．

	***  ARIMA Model Fitted to Series icecream ***

ここからレポートの内容を順に説明していく．

Model :  1 1 1 
Period:  1 

ModelはARIMAモデルの次数，階数を表しており， 左から順に自己回帰の次数，差分の階数，移動平均の次数になっている． このレポートの場合ARIMA(1,1,1)モデルを表している． Periodは季節成分の調整についてその周期を表しており， このレポートの場合季節成分が周期1であるため， 季節成分なし，つまり季節成分の調整を行っていないことを表している．

Coefficients:
	AR : -0.82439
	MA : -0.6889

CofficientsはARIMAモデルの係数を表している． ARは自己回帰の係数でMAは移動平均の係数である． このレポートの場合1次の自己回帰の係数は-0.82439， 1次の移動平均の係数は-0.6889であることを表している．

AIC: 2613.7006

AICは赤池情報量基準（Akaike's Information Criterion）という， モデルのよさの指標である．AICは低いほうがよいモデルであるとされる． このレポートの場合AICは2613.7006であることを表している．

4.4.2 残差のACFとPACFのグラフの読み方

残差とは，における実測値とモデルを用いて計算されたにおける推定値との差のことである． モデルの定義から残差の系列は定常であることが要求される． そこでACFのグラフをチェックする．

ACFの値が95%信頼区間を表す点線の内側に入っていれば（2,3個外側でも構わない），残差の系列が定常であるといえる． またACFのグラフがlagが増えるにしたがってあまり減少していかない場合には，残差にトレンドが残っている可能性がある．さらに一定の周期で増減を繰り返している場合には，残差に季節変動が残っている可能性がある．

先ほど立ち上がったグラフのウインドウにおいてPage1のタブをクリックすると，残差のACFとPACFのグラフが出力されている．

図x.ARIMA(1,1,1)モデルの残差のACFとPACF

ACFのグラフを見てみると，ACFの値が8つ95%信頼区間から外れている（lag0を除く）． これは，残差の系列が定常でないことを示している． また，明らかな1年周期の増減がみられる． これは，残差が季節変動していることを示している． この2点から，このモデルは再考の余地があることがわかる．

PACFからもモデルに再考の余地があることが読み取れる．

4.4.3 予測のグラフの読み方

先ほどのグラフのウインドウにおいてPage2のタブをクリックすると，予測のグラフが出力されている．

図x.ARIMA(1,1,1)モデルを用いた予測のグラフ

横軸は年（小数表示），縦軸が来客数の予測値である． また，実線が予測値，点線は予測の95%信頼区間である．

予測結果を見てみると，上昇のトレンドも季節変動も考慮されておらず，非常に疑わしい結果になっている． このことからもモデルを再考すべきという結論を得る．

4.5 ARIMAモデルの再構築

4.5.1 設定

この節では，レポートとグラフを活用して，さらによいモデルを選択する．

前節で"icecream"のモデルの検証をしたところ，残差に季節変動が残ってしまっていることがわかった． そこで，季節変動を除去したARIMAモデルを構築する．

図x.ARIMAモデルウィンドウ(Model)

モデルを構築する際に，季節変動の周期が1年であることを考慮して，[Seasonality]は[Monthly(12)]を選択する． 出力結果を上書きしないよう，別の名前"icecream.result2"として保存する．

図x.ARIMAモデルウィンドウ(Diagnostics)

図x.ARIMAモデルウィンドウ(Forecast)

出力結果を上書きしないよう，それぞれ別の名前"icecream.diagnostics2"，"icecream.forecast2"として保存する． 設定が終わったら[OK]をクリックする．

4.5.2 モデルの評価

レポートとグラフを用いて，再構築したモデルを評価してみる．

まず，レポートの核の部分を抜粋する．

	***  ARIMA Model Fitted to Series icecream ***

Model :  1 1 1 
Period:  12

Coefficients:
	AR : 0.01038
	MA : -0.04794

AIC: 2262.43397

ModelはARIMA(1,1,1)で先ほどと変わらないが， Periodが12になっている．これは年周期の季節成分の調整を行っていることを表している． Cofficientsから1次の自己回帰の係数は0.01038，1次の移動平均の係数は-0.04794に変わったことを表している． AICの値は2262.43397になり，先ほどのモデルよりも値が小さかった．これは，先ほどのモデルよりも当てはまりがよいことを示している．

次に，グラフを見てみる．

図x.季節変動を考慮したARIMA(1,1,1)モデルの残差のACFとPACF

残差のACFを見ると，95%信頼区間の外側になったのは2つだけであった．これは，残差の系列が定常であることを示している． また，顕著な季節変動もなくなった．ACFのグラフからも先ほどのモデルよりも当てはまりがよいことがわかる．

図x.季節変動を考慮したARIMA(1,1,1)モデルを用いた予測のグラフ

予測結果を見てみると，季節変動が考慮されている様子がわかる．

以上の3点から，モデルの当てはまりは季節変動を考慮したことによってよくなったことがわかる． もっともよいモデルを探すには，以上の観点を押さえたうえで， 次数や季節変動の有無を変えてみる探索的な方法で行う（コマンドラインを使えば自動でパラメータの同定を行うこともできる）． ちなみに，下の図は季節変動を考慮したARIMA(2,2,2)モデル（AICは1836.55213とさらにモデルの当てはまりがよくなっている）を用いて5年先まで予測した結果である． 季節変動に加えて増加のトレンドも反映されている．

図x.季節変動を考慮したARIMA(2,2,2)モデルを用いた予測のグラフ

"icecream"はデータ数が少ないため，ある程度次数を上げるとエラーが起こる（季節変動を考慮すると情報を多く失う）． そのため，エラーが起こるまで次数を上げれば上げるほどAICが下がる． しかし，一般的に次数を上げるとあるところで，AICは減少から増加に転じる．これはAICがパラメータの増加に対してペナルティーを課しているためである． 実際に"co2"を使ってAICの推移を調べたのが次の表である．なお，すべて季節変動を考慮している．

以上のような手順で，ARIMAモデルを用いた時系列分析を行っていく．