ボックス・ジェンキンズ法（１）

最初の時系列データ

前の授業では自己回帰という見方を紹介したので、今回はその自己回帰という観点に立てば将来予測ができるという方向へ話を発展させていこう。

最初に使用するデータファイル“IntroData-2015-M4.RData”を先に読み込んでおく（パワーポイントを参照のこと）。これはRの途中結果を保存したRDataファイルなので、ただ読み込むだけでよい–ダブルクリックするとRが自動的に起動され、作業用ディレクトリーもファイル保存場所に設定される。すでに保存されているデータや計算結果なども（丸ごと）復元される。

ただ、自動的にRが起動されるのは便利であるが、作業はRStudioで行いたい。この場合は、一度RStudioを立ち上げてから、RDataファイルを読む。メニューから“Session：Load Workspace”を選ぶ。あるいは、下のようにloadコマンドを使う。

load(file="IntroData-2015-M4.RData") # こちらを実行する

（注） （エクセルで作成した）CSVファイルは読み込んだ内容をRの中で記憶しなければならない（当たり前！）。故に、次のように（たとえば）データ全体に“data0”という名前をつける。

data0 <- read.csv(file="IntroData-2015-M4.csv") #これは実行しない

さらに、CSVファイルならRで読み込む前にどんな個別データがファイルに保存されているのかを、（ファイルを開いて）確認することも大事になる。

他方、RDataファイルの中身は全てRの作業結果である。それ故、どんな名前のデータが記憶されていたかをまず確認する。

ls()

## [1] "sdp"  "varp" "x1"   "x2"   "x3"   "x4"

そうすると、“x1”から順に“x4”までのデータがいま読み込まれたと分かる。

ここでデータ名を入力すればそれぞれの値が画面に表示されるのだが、ウッカリやると、たとえば100万個のデータが画面に流れ始める。最初の何個かをみるのが定石である。

（注） ウッカリ、巨大なデータを表示し始めたらescキーを押せばよい。

head(x1); head(x2);head(x3);head(x4)

## [1] 1.331936 1.705332 1.026151 2.717362 4.774558 7.698804

## [1] 0.7430316 1.8510006 3.0250518 5.2897208 5.4658500 6.1639571

## [1] 100.00000  98.03010  96.50879  97.47062  99.74165 101.35966

## [1]  0.0000000  0.4334259 -0.2104017 -1.0122456 -3.5235661 -5.1366678

どれも1種類ずつの系列である。とすれば、データのサイズを確かめるには長さ“length”を使う。縦横の大きさを“dim”で確認しようとしても“NULL”となり「出ません」という結果になる。

length(x1); length(x2); length(x3); length(x4)

## [1] 50

## [1] 50

## [1] 51

## [1] 51

ちなみにセミコロン記号で区切れば、複数のコマンドを同じ行に入力して一度に実行させることができる。

いずれの系列も50個前後で小さなデータである。

たとえば“x1”を図に描いてみよう。

plot(x1)

このように折れ線グラフが出てくる。

Rで描くグラフの標準形式は散布図である — 一変量の場合は横軸をデータ番号にとった散布図である。折れ線グラフが描かれるということは、データが既に時系列データになっている証拠である。

日付がつけられてなくとも時系列データになっているかどうかを確かめるには、次のようにするデータ型を確認する（注：Rには多くのパッケージがあり、それに応じて多くのデータ型がある。あまり意識する必要はない）。

class(x1)

## [1] "ts"

x1

## Time Series:
## Start = 1 
## End = 50 
## Frequency = 1 
##  [1]  1.331936  1.705332  1.026151  2.717362  4.774558  7.698804  6.211663
##  [8]  8.328886  9.045799  9.846150 11.221712 12.054022 12.032046 13.865109
## [15] 14.923664 17.608844 15.446357 17.631703 20.030811 19.208963 20.119859
## [22] 22.506899 23.718710 24.730306 25.611093 26.885922 25.959668 25.629842
## [29] 28.717924 28.789899 30.778027 31.412460 33.544502 35.246418 33.742872
## [36] 35.642588 34.480994 38.263227 37.470402 39.273885 39.394195 41.198633
## [43] 43.704990 42.850334 46.238908 44.649614 47.243412 47.431046 49.324990
## [50] 51.169319

Rに相当慣れてくるまでは「データ型」にまで目が向かないものである。上のコマンドで得られる“ts”という結果は、データ系列“x1”は時系列データ(ts)であるという意味である。但し、時点番号はつけられているものの、カレンダーの中の日付はない。

データ系列“x1”には明らかな増加傾向がある。将来に向けて「延ばしていく」という単純な予測も自然であるように感じる。が、「延ばしていく」という予測は、「これまで同じだったから、これからも同じであろう」という見方に比べれば、予測のロジックのうえで説得力にかける。

この疑問は、毎期に見られるデータ系列の変化を前期差をとることで氷解される。

plot(diff(x1))

mean(diff(x1))

## [1] 1.017089

sd(diff(x1))

## [1] 1.448685

前期差のグラフをみると、データ“x1”は平均1.0で増えており、増え方の標準偏差は1.45、つまりデータ“x1”は厳密に1ずつ増えてきたわけではないが、概ね1くらいずつ増えてきたとは言えるのであり、増え方は大体$1 \pm 1.45$の範囲に収まる。故に、これから将来にかけても（これまでどおり）毎期1程度の増え方で増えていくであろう。そう予測する。これが「延ばしていく」という発想のロジックである。

統計分析でいう「将来予測」は、統計的に安定して観察されてきた特徴を、将来に向けて伸ばしていく作業である。

授業ノートに登場するデータ例

前節でみた時系列データの中には毎期平均して1ずつ増えてきたデータが二つあった。しかし、その増え方、つまり変化のパターンを前期差のグラフでみると、どうも違うのではないか。そういう気づき方をした人はいなかっただろうか。これは将来予測のうえで非常に重要な点である。

今度は授業ノートにも出てくるデータを見ながら、変化のパターンを掘り下げて調べることにする。データファイルはやはりRDataであるからloadする。

load(file="TextData-2015.RData")
ls()

##  [1] "Data1" "Data2" "Data3" "Data4" "sdp"   "varp"  "x1"    "x2"   
##  [9] "x3"    "x4"

授業ノートに登場するデータはすべてRDataファイル“TextData-2015.RData”に保存されている。これをloadコマンドで読み込んで、データ一覧す。そうすると、先ほど読みこんだ“x1”から“x4”までに系列に加えて“Data1”から“Data4”までが新たに加わっている。これがいま読み込んだデータである。

plot(diff(Data1))

mean(diff(Data1))

## [1] -0.04717719

sd(diff(Data1))

## [1] 0.9715799

dx1 <- diff(Data1)
lag.plot(dx1)

lag.plot(dx1,do.lines=F)

lag.plot(dx1,do.lines=F,diag=F)

lag.plot(dx1,lags=4,do.lines=F,diag=F)

acf(dx1,lag.max=20)

自己相関図・偏自己相関図

自己相関図については頭に入っただろうか？少し以前までは、まず変化のパターンを判別し、どんな型に分類されるかを図で読み取ってから、当てはめる予測モデルを決め、それを推定していた。かなり煩雑であり、相当の分析経験を必要としたのだが、最近はソフトウェア自体が自動的に予測モデルを探索する機能を備えつつある — Rのパッケージ“forecast”もそうである。

とはいえ、まったくの「機械任せ」で予測の基本ツールを知らないと、レポート作成・説明において自信が持てず、説明が不安定になる。最小限の知識・理解はあったほうがよい。

そのため、データファイルを変えて実習することにしよう。一度、ここで記憶したデータ系列をすべて消去しておく。Windowsでは＜その他→全てのオブジェクトの消去＞を選ぶ。RStudioを使用している場合は、＜Session→Clear Workspace＞をメニューから選ぶ。

なお、今回の実習ではパッケージ“forecast”は使わずともよい。最初から搭載されている機能だけで実習を進めることにしよう。

データファイル”FiveSeries.RData”を読み込む。以下のようにしてもよいし、データファイルをダブルクリックするとRのコマンド画面が立ち上がる。RStudioでは＜Session → Load Workspace＞を選ぶ。

load(file="FiveSeries.RData")
ls()

##  [1] "data.ima" "data1"    "Data1"    "data2"    "Data2"    "data3"   
##  [7] "Data3"    "data4"    "Data4"    "data5"    "dx1"      "sdp"     
## [13] "varp"     "x1"       "x2"       "x3"       "x4"

データ系列”data.ima”の変化の特徴を調べてから、将来6期予測をする。但し、今回は1.同定、2．推定、3．診断、4．予測というボックス・ジェンキンズ法を構成する四つのステップのうち診断ステップは省略する。

plot(data.ima)

データ系列“data.ima”は増加傾向と減少傾向が交代で現れている。このように、期間全体で一貫した傾向ではないが、相当期間にわたって交代で現れる傾向を『確率的なトレンドがある』という言い方で表現している。 確定的であれ、確率的であれ、トレンド要素があるとき、データの自己相関をみても変化のパターンは把握できない。

まず“data.ima”について、4期ラグまでの散布図を描いてみよう。

lag.plot(data.ima,lags=4,do.lines=F,diag=F)

すべてプラスの相関がある。元のデータには増加局面、減少局面の両方があるにもかかわらずである。

自己相関図にしても同じである。

acf(data.ima)

自己相関はプラスからゆっくりと低下し、やがてマイナスになり、ゆっくりと下がっている。

このような自己相関図は、高止まり型と呼んでおり、データの前期差をとって変動パターンを見るべきことを教えている。

acf(diff(data.ima))

自己相関図をコレログラムと呼ぶ専門家も多いが、今回得られた図をみるとラグ1の相関が高い — ラグ0は常に1である。前後に比べて突出した箇所をスパイクというが、図で見られるスパイクは1本だけである。

このような自己相関図をカット型と呼んでいる。実は、自己相関図がカット型であるデータは、前回授業でとりあげた自己回帰分析では正確に予測できない。これが今までに得られている理論的帰結である — 指数平滑法が当てはまるのは今回のような変動パターンを示すデータに限られる。

更に色々な変動パターンをみていこう。

データ“data1”はどうだろう。

plot(data1)

plot(diff(data1))

mean(diff(data1)); sd(diff(data1))

## [1] 0.305121

## [1] 1.28738

“data1”は全体として増加している。ということは、前期差の平均はプラスのはずである。

データ“data1”には明らかなトレンドが含まれるので、元データの自己相関図を描いても高止まり型になると予想できる。

【クイズ】

トレンドが混じっていると、なぜこんな結果になるのだろうか？

データ“data1”の前期差をとってから自己相関図を描く。

acf(diff(data1))

pacf(diff(data1))

“data1”の前期差の自己相関図は、どこかに突出して高い相関があるカット型ではない。むしろ波の形を描きながらユックリと小さくなっている。一度小さくなっても、また高い相関が現れる。そんな形である。これを減衰型と呼んでいる。

二番目のコマンドは“acf”と似ているが、これは偏自己相関図を描くコマンドである。描かれた図の読み方は自己相関図と同じであるが、偏自己相関図では左端がラグ0ではなく、ラグ1が左端になっている。

理論の詳細はともかく、偏自己相関図が伝える意味合いは上につきる。

得られた偏自己相関図はラグ1とその右のラグ2にスパイクがある。他は青い線の範囲におさまり有意でない。こんな場合、2次の自己回帰モデルがデータに当てはまることということを図が伝えている（パワーポイントを参照）。

同定から推定・予測へ

データ“data1”の前期差については、自己相関図が減衰型、偏自己相関図が二本スパイクのカット型だった。こんな場合は、前期差について2次の自己回帰モデルが当てはまる。これを式で書けば下のようになる。

$$\Delta X_t -m = a_1 (\Delta X_{t-1}-m) + a_2 (\Delta X_{t-2}-m)$$

前期差である$\Delta X_t$に対して自己回帰式を当てはめている。但し、$m$は前期差の平均値であって、これがプラスなら元の値“data1”は全体としては増加している、マイナスなら減少していることになる。前期差の平均がゼロなら、元の値は期間を通しておおむね横ばいであるはずだ。データ“data1”は期間を通して増加しているので、前期差の平均値“m”にはプラスの値を想定するケースである。

まず、前期差についての自己回帰式を推定してみよう。パッケージ“forecast”でなく、標準でRに搭載されているコマンドを使ってみる—今後はあまり使わないが。

kekka1 <- arima(data1,order=c(2,1,0))

推定の計算法は回帰分析用いる単純な最小二乗法ではない。

予測式が推定されたあと、推定結果に基づいた予測計算には“predict”コマンドを用いる。

yosoku1<-predict(kekka1,n.ahead=6)
ts.plot(data1,yosoku1$pred,lty=c(1,2),col=c(1,2))

データ“data1”の最終期は101だから、その先6期間予測は102期から107期までとなる。予測値は期待値（=中位予測）と予測値の標準誤差が計算される。

グラフに示されている予測経路をみると、概ね横ばいであり、全体に見られる増加傾向は繁栄されていない。実は、何も指定せずにarimaコマンドで推定すると前期差の平均はゼロであると前提される — この辺をどう修正すれば予測らしい予測になるのか。これも次回に議論することにしよう。

以上、データ“data1”を例にとって予測計算まで行ってみたが、データに含まれる変動パターンから、どんな予測モデルがデータに当てはまっているかを判別するステップが出発点だった。判別された予測モデルをarimaコマンドの中のorderに指定するのだが、ここがボックス・ジェンキンズ法という予測実務の勘所になる。

kekka2 <- arima(data.ima,order=c(0,1,1))
yosoku2 <- predict(kekka2,n.ahead=6)
ts.plot(data.ima,yosoku2$pred,lty=c(1,2),col=c(1,2))

data1.sub <- window(data1,end=95)

kekka1.sub <- arima(data1.sub,order=c(2,1,0))
yosoku1.sub <- predict(kekka1.sub,n.ahead=6)

ts.plot(data1.sub,yosoku1.sub$pred,col=c(1,2))
ts.plot(data1,yosoku1.sub$pred,col=c(1,2))