時系列の基本成分（TCSI）

住宅着工・月次データを見る＜テキストの再確認＞

　データ入力

Rのコマンドでデータファイルを読み込むときのルーティンは以下のとおりである。

まず（必ずしも不可欠というわけではないが）作業用ディレクトリーを所定のフォルダーに移した後、以下の順番で実行する。

ファイル一覧でファイル名を確認
データファイルを読みこむ。
読み込んだデータの大きさ（行数×列数）を確認
最初の何行かを見る

まず住宅着工戸数を保存したデータファイルを読み込む。データはデータベースで管理されていることが多いが、分析する際はファイルの中身をエクセルで表示でき、すべての統計分析ソフトで読み込めるCSVファイルで保存することが多い（注：　CSVとは“Comma Separated Values”の略称で、エディターでも開くことができる）。ファイルの第１行は個別データの名前(header)であるのが普通である。

dir()

## [1] "AustralianBeerProduction.xls"    "forecast-m2-2015-ansi.html"     
## [3] "forecast-m2-2015-ansi.rmd"       "forecast-m2-2015-utf8.html"     
## [5] "forecast-m2-2015-utf8.rmd"       "HousingStart-196501-201207.csv" 
## [7] "hs-198501-201408.csv"            "TrendRegressionCons.RData"      
## [9] "住宅着工戸数-196501-201408.xlsx"

data0 <- read.csv(file="HousingStart-196501-201207.csv",header=T)
dim(data0)

## [1] 571   3

head(data0)

##   year month    hs
## 1 1965     1 54772
## 2 1965     2 62092
## 3 1965     3 61206
## 4 1965     4 72861
## 5 1965     5 64390
## 6 1965     6 73223

【補足】もし数行ではなくもっと表示させたい場合、また最後の数行をみたい場合は次のようにする。

head(data0,15)

##    year month    hs
## 1  1965     1 54772
## 2  1965     2 62092
## 3  1965     3 61206
## 4  1965     4 72861
## 5  1965     5 64390
## 6  1965     6 73223
## 7  1965     7 74464
## 8  1965     8 69585
## 9  1965     9 72798
## 10 1965    10 84160
## 11 1965    11 83541
## 12 1965    12 69504
## 13 1966     1 54260
## 14 1966     2 62383
## 15 1966     3 63939

tail(data0,3)

##     year month    hs
## 569 2012     5 69638
## 570 2012     6 72566
## 571 2012     7 75421

住宅着工データを時系列データ(ts)にする

データファイルには三つの個別データが保存されていた。データセット全体にはdata0という名前がいま付けられたが、分析対象は3列目の住宅着工戸数である。しかし、その3列目をとって図に描いてみても、時系列として見やすいグラフにはならない。

時系列データとしての扱いを便利にするために3列目をとりだし、日付をつけ、時系列(ts)型の変数として別データとしてもっておくことにする。

hs <- data0[,3]
plot(hs,main="住宅着工戸数")

hs <- ts(hs,freq=12,start=c(1965,1))
plot(hs)

時系列データを図に描くと、自動的に折れ線グラフが選ばれ、横軸には日付が表示される。

この住宅着工データは月次データであるので、季節性があるはずである。ところがデータ期間が長いためグラフでは明瞭にデータの特徴が現れていない。そこで、1985年1月から1990年12月まで、更に2000年1月から2007年12月までを切り取って、別に名前をつけて保存しよう。

時系列データの一部期間を切り取るにはwindowコマンドを使う。

hs.sub1 <- window(hs,start=c(2000,1),end=c(2007,12))
hs.sub2 <- window(hs,start=c(1985,1),end=c(1990,12))
plot(hs.sub2)

length(hs.sub2)

## [1] 72

期間を限定して観察することで、住宅着工戸数は毎年同じパターンで変動を反復していることがみてとれるだろう。これが季節性である。経済・経営領域の月次系列、四半期系列では、大方の場合、季節性が混じっているものである。

直線傾向線（トレンド）の当てはめ

住宅着工戸数の月次系列には季節性の他に、景気変動に伴う循環、日本経済の潜在成長力を反映する長期トレンドといった成分が混在している。

時系列データは四つの基本成分、即ち傾向(T)成分、循環(C)成分、季節(S)成分、不規則(I)成分から構成されていると見るのが基本的なアプローチだが、与えられた時系列データにどんなトレンドが含まれているかを検証するには幾つかの方法がある。

時系列データからトレンド(T)成分を、あるいはトレンド・循環(TC)成分のみをとり出したり、また季節成分を除いてTCI成分にする季節調整を、通常、フィルター、もしくは（ニュアンスはやや違うが）平滑化と呼んでいる。

最も簡便な方法は、特定の型のトレンドを回帰分析で当てはめる方法である。たとえば直線傾向線を当てはめるには以下のようにする。対象としては上で保存した部分系列“hs.sub2”にしよう。データの長さは72か月である。

time <- seq(1,72)
kekka1 <- lm(hs.sub2 ~ time)
yosoku1 <- predict(kekka1) # 結果から計算値をとりだす
head(yosoku1) # 日付がついていない

##        1        2        3        4        5        6 
## 104377.6 105087.3 105797.1 106506.9 107216.6 107926.4

yosoku1 <- ts(yosoku1,freq=12,start=c(1985,1))
plot(hs.sub2)
lines(yosoku1,col="red",lty=2,lwd=2)

Rで回帰分析を行ったあと、元のデータを散布図に描き、その図の中に回帰直線を描き足すのは簡単だった—abline(回帰分析の結果名)でいい。

ところが時系列データの場合、計算値に日付けがついていないので、上のように計算値に日付をつける必要がある。

図の直線傾向線は、住宅着工戸数の中のT成分である。ということは、データの値からT成分を引いた値、つまり回帰分析の残差が残りのCSI成分だということになる。回帰分析の残差の平均は常にゼロである点にも注意せよ。簡単に言えば、住宅着工戸数の平均的な水準をトレンドとしたあと、数年期間のプラスマイナスの変動をC成分、1年期間のプラスマイナスをS成分、毎月のランダムなプラスマイナスをI成分として見るわけである。

zansa1 <- residuals(kekka1)　# 結果から残差をとりだす
zansa1 <- ts(zansa1,f=12,s=c(1985,1))　#　残差に日付け
par(mfrow=c(2,1))　#　2行、1列の段グラフを描く
plot(yosoku1)
plot(zansa1)
abline(h=0)　#　グラフに水平線を描き足す（目盛ゼロ）

sum(zansa1)　#　確かに残差の合計（平均）はゼロか

## [1] -2.467004e-10

移動平均で季節成分を消す

直線トレンドを回帰分析で当てはめる方法は、簡単ではあるが、トレンドの型を最初から直線に前提してもよいかという問題が残る。あくまでも簡便法として使うべきだろう。

これに対して（計算自体はシンプルだが）移動平均法は、季節調整の基礎にもなっており、広く活用されている平滑法である。

移動平均とは、定められた数のデータの平均値を日付に対応させていく方法である。平均計算は単純平均であったり、加重平均であったりする。

それでは、系列“hs.sub2”の3か月移動平均を求めよう。移動平均は“filter”コマンドで行う。

ma3.sub2 <- filter(hs.sub2,rep(1/3,3))
head(hs.sub2)

## [1]  77527  90942 103955 116794  99660 112363

head(ma3.sub2)

## [1]       NA  90808.0 103897.0 106803.0 109605.7 106377.7

ts.plot(hs.sub2,ma3.sub2,lty=c(1,2))

図の実線が元のデータ、点線が３か月移動平均値である。

filterコマンドの括弧内にあるrep(1/3,3)は、分数\(1/3\)を3個並べた列（＝Rでいうベクトル）を表し、コマンドの意味合いは時系列データを\(X_t\)として、以下の平均値を求めているわけである。

\[ MA_t = \frac{1}{3}\times X_{t-1} + \frac{1}{3}\times X_{t} + \frac{1}{3}\times X_{t+1} \]

元データの最初の3項の平均値が移動平均値の1番目の値と一致していることを確かめよ。

上は中央移動平均である。つまり3か月分を平均した平均値をとった3月の中央に対応付ける。

これとは別に後方移動平均もある。これは、例えば株価の25日移動平均や90日移動平均線という形で利用されている。

ma3.sub2.kouhou <- filter(hs.sub2,rep(1/3,3),sides=1)
head(hs.sub2)

## [1]  77527  90942 103955 116794  99660 112363

head(ma3.sub2.kouhou)

## [1]       NA       NA  90808.0 103897.0 106803.0 109605.7

3か月分のデータをとって最後の月に対応付けるので、最初の2月分の平均値はNA、つまり欠落項になる—中央移動平均では1月分の欠落項がでる。

3項移動平均は、特に不規則成分の除去を目的に、しばしば利用されている。

では、住宅着工戸数のデータから季節成分をとりさることにしよう。それには12カ月移動平均をとればよい。しかしながら、偶数項をとった中央移動平均をとると、平均値の対応月が半月ずれてしまう。たとえば偶数である４月移動平均をとると、1月から4月までの平均値が2.5月という時点に対応してしまう。この時点のずれを修正するため2項移動平均を行っておく必要がある。即ち、2段階の移動平均が必要になるのであるが、この計算を1回で行うこともできる。それには、以下のような加重移動平均を実行すればよい。

wgt13 <- c(1/24,rep(1/12,11),1/24)　#　13項加重移動平均のウェイト
ma13.sub2 <- filter(hs.sub2,filter=wgt13,side=2)
ts.plot(hs.sub2,ma13.sub2,type="l",lty=c(1,2))

式で書くと以下のようになる。元のデータを\(X_t\)、現時点を\(t\)時点とする。

\[ MA_t = \frac{1}{24}\times X_{t-6} + \frac{1}{12} \times X_{t-5} + \cdots + \frac{1}{12}\times X_{t+5} + \frac{1}{24}\times X_{t+6} \]

図をみると、13項移動平均によって季節成分が消え去り、平滑化されていることが分かるだろう。なお、13項の移動平均をとっているから、元のデータに含まれていたはずの不規則成分も相互に相殺されているはずである。故に、13項移動平均値は元のデータにあったTC成分、つまり傾向＋循環成分であると解釈できる。

　C、S、I成分への簡単な分解

月次系列や四半期系列が与えられて、明らかに季節成分が含まれている場合、まずは季節成分を除去して平滑化し、データの中の基調的な変化を見ることが時系列データ分析の第一歩である。それにはfilterコマンドが有用なのだが、やや手数がかかって面倒である。

コマンド“decompose”は上で述べた計算を一度で実行するので、普通はこちらを使えばいいだろう。

bunkai <- decompose(hs.sub2)
plot(bunkai)

# どんな結果が計算されたのか。諸々の結果の名前を確認する
names(bunkai)

## [1] "x"        "seasonal" "trend"    "random"   "figure"   "type"

移動平均を行うと必ずデータの両端部分に欠落項が生まれる。

ただ、季節成分は同じパターンが固定的に反復されるので、分解された図でも全期間にわたって描画されている。しかし、データの始期、終期の部分は、本来は欠落項になっているので、TC成分とI成分が分離されていない。そのため図の両端は欠けている。これは移動平均という方法の欠点であるが、最近では新たなフィルタリングが開発・提案されており、足元の変化がわからないという状況ではない。

授業では、次回以降、自己回帰分析へと方向を転じて、ボックス・ジェンキンズによるARIMA分析へ進んでいくので、発展した移動平均法はここまでにする。

【参考】ホルトウィンタース法による予測法　※田中『Rによる時系列分析入門』（pp111～121）

hw.sub2 <- HoltWinters(hs.sub2,seasonal="mult")
fitted.sub2 <- hw.sub2$fitted[,1]
ts.plot(hs.sub2,fitted.sub2,type="l",lty=c(1,2))
yosoku.sub2 <- predict(hw.sub2,n.ahead=12)
ts.plot(hs.sub2,yosoku.sub2,type="l",lty=c(1,2))

　練習問題

以下の作業を行ってください。

上で作成した住宅着工データのサンプル“hs.sub1”について、同じ作業を行ってください。
TC成分は将来予測できそうですか？あなたならどうしますか？季節成分はサンプル期間の外に伸ばせそうですか？不規則成分についてはどう思いますか？
ファイル“TrendRegressionCons.RData”を使って、授業ノート16ページの練習問題に取り組んでください。
ファイル“AustralianBeerProduction.xls”を使って、オーストラリアの月次ビール生産量の変化について分析してください。

【ファイルを開くときの補足】

ファイル名の修飾子がRDataであるファイルは、Rで保存されているデータ内容、計算結果等をすべて書き込んだファイルである。RDataファイルをダブルクリックすると自動的にRが起動される。どんな変数が保存されているかを確かめるにはlsコマンドを使う。この点は改めて後述する。
エクセルで作成したxlsファイルを読むには、CSVファイルで保存しなおすのが簡単。パッケージ“xlsxを利用すると、エクセルで作成したワークブックファイル（*.xlsx）から読み込むワークシートを指定して読み込めるが、ここでは省略する。