時系列解析とは？時間の多項式トレンド(polynomial trends)を使った基本姿勢、計算方法とRコードの解説

そういえば時系列に関する話、あまり書いていないなと思ったので記事にしていこうと思います。

時系列とは？

まず、簡単に時系列とは何かっていう話ですが、例えば株価とかを考えてもらうと、時間の経過とともに何か値が連動して変化していますよね。
f:id:bananarian:20181011013926p:plain

こういうのを時系列データと呼びます。

時系列解析の種類

時系列解析では大きく分けて二つの動きを考えます。

トレンド

「上がったり下がったりすれど、大きく見ると時間と共にじわじわ増えているよね～」だとか、「じわじわ下がっているね～」だとか考える時、人はトレンドを考えています。

要は何か一定の周期があったりするわけではないけども、大雑把に見ると時間と共に流れる推移が見えるような場合ってありますよね。その傾向は大抵そのまま続くものです。

このような大きく見た時の傾向をトレンドと呼びます。

周期

一方、何かしらの影響が背後にあって、一定の範囲で周期性があるような動きを周期と呼びます。
例えば、日本の気温の推移を考えてもらえばいいと思うんですけど、

春から夏にかけてどんどん気温が上がっていき、秋、冬にかけて気温が下がっていく、そしてまた春になって同じような動きをする

こんな動きを周期と呼んだりします。

時系列解析ではこの周期とトレンドの二つの観点からものを考えていきます。
例えば気象庁のサイトには日本の平均気温のグラフとして下のようなものがあったりします。

f:id:bananarian:20181112113246p:plain
(気象庁 | 日本の年平均気温より)

この赤い線がトレンドで、上昇傾向にあります。

時系列の基本スタンス

確かに直感的にも時間を伴うデータは大きなトレンドと周期に分けられそうですよね(?)。

とよく説明されるわけですが、正直思いとしては

せめて周期(適当な規則性)がないと、予測も何も出来ないからそれくらいの仮定は設けさせてくれだと思ってます。

だってトレンド(大きな傾向)もなくて、周期(適当な規則)もない各時間ごとに1回こっきりしか観測されないデータを予測するなんてそもそも無理ですよね。だから無理のないやつを予測対象におきたいわけです。

そんな思いがあるわけなので、時系列解析を行う場合には大原則として注意点があって、

「そもそもこのデータには、大きな傾向や規則が生じうるのだろうか」

ということを考える必要があります。

基本のトレンドモデル

前置きが長くなりましたが、今回はそんな時系列解析の内、トレンドを表現するモデルを説明していこうと思います。

トレンドの予測にも色々あって、例えば企業の利益や株価等の経済データを考えてもらうと、これらのデータは大抵は上昇トレンドを持っています。これは何故かというと、基本的に経済は時間と共に成長を続けているはずなので、その成長に連動して、各企業の利益や株価もじわじわ増えていくことが関係しています。

※当然例外はあります

このように背後の要因がわかっているような場合は、経済成長を示すような指標を使うことで株価のトレンドなどを推定することが可能です。

しかし、当然何かトレンドはありそうだけど、背後の要因が何か検討もつかないなんてこともあり得ますよね。

でも私たちは時系列解析を行う以上、確実にトレンドに関係しそうな情報を１つ持っているわけです。それが

時間

です。

そこで、(本当はもっと本質的に関係しそうな情報を使うのがベストだが)データに関して何の情報も持ち合わせていないときに、とりあえず時間の推移を軸に、多項式近似モデルを作って、トレンドを推定してやろうというトレンドモデルが存在します。

それが次のようなモデルです。

時点 $t$ の予測したいデータを $y_t$ とおく。その時パラメータ $\alpha_0,....\alpha_q$ を使って次のような $q$ 次多項式でトレンドを近似する。

$y_t=\alpha_0+\alpha_1 t+\alpha_2 t^2+...+\alpha_q t^q$

当然、完全にこの式にピッタリフィットするはずがないので、確率的な誤差 $\epsilon_t$ も与えておきます。

$y_t=\alpha_0+\alpha_1 t+\alpha_2 t^2+...+\alpha_q t^q+\epsilon_t$

ちなみに、先ほどの気象庁のグラフのトレンドはこれの $q=1$ を利用しています。1次関数が右上がりになっていましたよね。

簡単そう

なんだ！よく見たらただの線形回帰じゃないか！簡単そうですね！！と思うかもしれませんが、確かに $q=1$ なら線形回帰ですが、 $q≧2$ ではそう簡単にいかないことに注意してください。

線形回帰分析の仮定として説明変数間の無相関仮定があったはずです。しかし、 $t$ と $t^3$ ってどう考えても相関ありますよね。

そこで、先ほどの式を無相関な説明変数が並ぶ形に書き換える必要が生じます。

無相関な関係

そこで、次のような条件を満たす $\phi_{iT}(t)$ に説明変数を書き換えるという操作を行います。

$\sum_{t=1}^{T}\phi_{iT}(t)\phi_{kT}(t)=0$
ただし $i≠k$ , $i,k=0,1,....T-1$

これは直交条件と呼ばれていて、要は説明変数間はこの仮定を満たす必要があると思ってください。

$\phi_{kT}(t)$ を $t$ を使って表すために、次の多項式で表現することにします。
$\phi_{kT}(t)=t^k+X_{k-1}t^{k-1}+....+C_1 t+C_0$

この $C_?$ を求めることが出来れば、説明変数の書き換えを行えますね。

ここで、直交条件が成り立つということは、次の条件が成り立つことがわかりますか？
$\sum_{t=1}^{T}\phi_{kT}(t)t^i=0$
ただし $i≠k$ , $i=0,1,....k-1$

要は先ほどの式に放り込んでみると、 $k$ 以外の $t$ の何とか乗を掛けて全部足せば0になるという形にならねばならないことがわかります。

よって、先ほどの関係式
$\phi_{kT}(t)=t^k+X_{k-1}t^{k-1}+....+C_1 t+C_0$

に対して両辺に $\sum_{t=1}^T t^i$ を掛けてやると
ただし $i≠k$ , $i=0,1,....k-1$

次のようになりますよね。

$C_0 \sum_{t=1}^T t^i+C_1 \sum_{t=1}^T t^{i+1}+....+C_{k-1} \sum_{t=1}^T t^{i+k-1}=-\sum_{t=1}^T t^{i+k}$
ただし $i≠k$ , $i=0,1,....k-1$

何かゴツゴツしていますが、まず $k=1$ とおきます。すると $i=0$ となるので式は次のようになります。

$TC_0=-\sum_{t=1}^T t$

よって $C_0=-\frac{T+1}{2}$ だとわかりますよね。

次に $k=2$ を考えると $C_0$ がわかっているので $C_1$ が求まります。以下同様に繰り返せば全ての $C_?$ が求まりますよね。

このもとまった $C_?$ を活用して、

$\phi_{1T}(t),....\phi_{qT}(t)$ を用意し、

$y_t=\alpha_0+\alpha_1 t+\alpha_2 t^2+...+\alpha_q t^q=\gamma_0+\gamma_1 \phi_{1T}(t)+.....+\gamma_q \phi_{qT}(t)$

という問題に書き換えることが出来ました！というわけです。これなら線形回帰を行うことが出来ます。

Rで実装

このトレンドモデル、正直かなり古いモデルなので恐らくRのパッケージにはありません。そこで今の手順をそのままコードに直してみました。

###値の設定、tは時間,degreeは多項式の次数
###ここはデータに合わせて変更する
t=1:30
degree=3
###
#expl_matに変換した説明変数を行列形式で突っ込む(tと次数の分だけ増える)
expl_mat=matrix(,length(t),degree)
for(tt in t){
	
	k=1:degree
	phi=c()
	Coef=c()

	i=(k-1)
	for(kk in k){
		ii=i[kk]
		ti_mat=matrix(,length(t),kk+1)
		for(iii in ii:(ii+kk)){
			ti_mat[,which(ii:(ii+kk)==iii)]=t^iii
		}
		ti_seq=colSums(ti_mat)
		if(length(ti_seq)==2){
			Coef[1]=-(length(t)+1)/2
		}
		else{
			mti_seq=ti_seq[-((length(ti_seq)-1):length(ti_seq))]
			Coef[kk]=(-ti_seq[length(ti_seq)]-sum(Coef*mti_seq))/ti_seq[length(ti_seq)-1]
		}
		phi[kk]=sum(c(1,rev(Coef))*(tt^rev(seq(0,kk,1))))
	}
	expl_mat[tt,]=phi
}

力作！笑

さて、一応合っているか確認してみましょう。もし、変換が正しければ次数さえ真のデータとあっていれば変換前と変換後で値が完全一致するはずです。

#真のモデルが3次多項式のデータを作成
y=0.3+0.5*t+0.5*t^2+0.5*t^3

##間省略、先ほどのコードにdegree=3を突っ込みexpl_matを作る
##浮動小数点の関係で小数点以下までは一致しないので丸め込み(round)して同じ値かどうか確認

> round(expl_mat%*%coef(lm(y~expl_mat))[2:4]+coef(lm(y~expl_mat))[1])==round(y)
      [,1]
 [1,] TRUE
 [2,] TRUE
 [3,] TRUE
 [4,] TRUE
 [5,] TRUE
 [6,] TRUE
 [7,] TRUE
 [8,] TRUE
 [9,] TRUE
[10,] TRUE
[11,] TRUE
[12,] TRUE
[13,] TRUE
[14,] TRUE
[15,] TRUE
[16,] TRUE
[17,] TRUE
[18,] TRUE
[19,] TRUE
[20,] TRUE
[21,] TRUE
[22,] TRUE
[23,] TRUE
[24,] TRUE
[25,] TRUE
[26,] TRUE
[27,] TRUE
[28,] TRUE
[29,] TRUE
[30,] TRUE

というわけで、正しいコードを書くことが出来ました。
これで簡単に多項式トレンドを考えることが出来ますね。

多項式トレンドの使い方

まず、あまり次数を増やすとフィットはしますが、解釈が難しくなるので、次数は小さい方が望ましいです。
最も解釈しやすいモデルが1次式のため、気象庁ページも1次式を採用しています。

普通、次数を決める際は複数の次数で回帰を行ってみてその結果で決めます。その辺りは今回は省略しますが、仮説検定を行う方法や、AICを使う方法、stepwise法等色々あります。

ということで、時間だけを使った多項式トレンドモデルの説明をしました。次回また別の時系列について説明しようと思います。