時系列データでお絵かきする - データサイエンティスト(仮)

はじめに

最近、趣味の時間があまり取れなかったので更新が滞ってました（汗）。

周期性をもつ（だろう）時系列データをお遊びで可視化してみました。やってみるとなかなか面白かったので、備忘録として残しておくことにしました。

どうやって可視化するのか

ある時点 $t$ での観測値を $y_t$ とするような時系列データ $\{y_t, t = 1, \cdots, T \}$ を考えます。 $y_t$ と1期前の観測値 $y_{t-1}$ を成分にもつ2次元ベクトル
${ \displaystyle \begin{eqnarray} \xi _t= \begin{pmatrix} y_t \\ y_{t-1} \\ \end{pmatrix} \tag{1} \end{eqnarray} }$
を考え、これを $t$ を変えてプロットします。今回、1期前だけを考えましたが、並べるのは何期前のものでもよいですし、さらに過去の観測値を並べて高次元のベクトルを考えても良いです。とにかく現時点と過去の観測データを幾つかチョイスし、それらを成分として構成したベクトルをプロットする、ということが具体的にやることです。

可視化

サンプルデータ1 : 単一の周期性をもつノイズのない系列

単純な場合として、観測値が $y_t = \cos(t)$ となっている場合を考えます。

# ライブラリ
if(!require(data.table)){
  install.packages("data.table", dependencies = TRUE)
}
require(data.table)
if(!require(dplyr)){
  install.packages("dplyr", dependencies = TRUE)
}
require(dplyr)

# サンプルデータを100期分作成
t <- seq(1, 100, 1)
x1 <- cos(t)
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

観測値をobsとし、obs_lagはtを1期さかのぼったデータが格納されています。

観測値をプロットすると、たしかに周期性をもつような振る舞いのデータになっています。

# ggplot
if(!require(ggplot2)){
  install.packages("ggplot2", dependencies = TRUE)
}
require(ggplot2)

# 観測値をプロット
ggplot(data = df, 
       aes(x = time, y = obs)
) + 
  geom_line()

f:id:tekenuko:20161110225451p:plain

次に、obsとobs_lagの散布図を描いてみます。

# プロット
ggplot(data = df, 
       aes(x = obs, y = obs_lag)
       ) + 
  geom_point(color = "hotpink")

結果は以下になります。
f:id:tekenuko:20161110231730p:plain

傾いた楕円のような形に見えるものが生成されました。このような形になることは、 $\sin^2(t) + \cos^2(t) = 1$ を $y_t$ と $y_{t-1}$ で表すことで式でも理解できますが、ここでは省略します。言いたいことは、周期性のある時系列データに関して、時点 $t$ と過去の時点の観測値を並べたものを見ると、なにか構造が見えることがあるということです。

傾いた楕円のようなものが生成されるのは、もとの時系列が単一の周期を持っていたからです。この振る舞いは、比較する観測値を変えても変わりません。比較として、 $y_t$ と2期前の観測値 $y_{t-2}$ を成分とするベクトルをプロットした場合も考えます。中心や傾き方は変わりますが、傾いた楕円のようなものが生成されます。

# 2期前の観測値を並べる
df[, obs_lag2 := lag(obs, 2)]
# プロット
ggplot(data = df, 
            aes(x = obs, y = obs_lag2)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161110231928p:plain

サンプルデータ2 : 単一の周期性をもつノイズがある系列

さきほどはきれいな形が見れましたが、ノイズが付与されたデータだとどうでしょうか。今度は、 $\cos(t)$ にホワイトノイズ $N(0, 0.1)$ が付与されたデータを考えてみます。こちらの場合にobsとobs_lagの散布図を描くと

t <- seq(1, 100, 1)
x1 <- cos(t) + rnorm(100, 0, 0.1)
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

# プロット
ggplot(data = df, 
            aes(x = obs, y = obs_lag)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161110233410p:plain

となります。目を凝らすと楕円に見えなくもないです。このように、ノイズがある場合は、なんらかの構造にノイズがのったものが生成されます。

注意 : トレンドがある場合

トレンド（長期的な上昇/下降傾向）があるような時系列データの場合は、差分系列にしてからベクトル化します。以下のような上昇トレンドをもつようなデータを考えます。

t <- seq(1, 100, 1)
x1 <- cos(t) + t + rnorm(100, 0, 0.1) 
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

# プロット
g <- ggplot(data = df,
            aes(x = time, y = obs)
) +
  geom_line()
ggsave("fig_cau.png", g)

f:id:tekenuko:20161110234038p:plain

この場合、obsとobs_lagの関係は以下になります。

ggplot(data = df, 
            aes(x = obs, y = obs_lag)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161110234152p:plain

楕円とは似ても似つかないものが生成されました。これは、 $y_t$ が周期関数以外の部分で $t$ にあらわに依存しているためです。そのため、 $\sin^2(t) + \cos^2(t) = 1$ を $y_t$ と $y_{t-1}$ で表すと、 $t$ の寄与が残ってしまい、 $y_t$ と $y_{t-1}$ の関係式（ $(y_t, y_{t-1})$ 平面での軌道に対応）が $t$ に依存する形になってしまうのです。よって、軌道が閉じなくなり、上のような直線っぽい振る舞いになります。

周期性そのものを可視化したい場合は、このトレンドは調整した後にベクトル化する必要があります。差分系列を取った後、ベクトル化したものを可視化すると以下のようになります。

# 差分系列を作成
df[, diff_obs := ifelse(t == 1, 0, obs - obs_lag)]
df[, diff_obs_lag := lag(diff_obs)]

# プロット
ggplot(data = df, 
            aes(x = diff_obs, y = diff_obs_lag)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161111000219p:plain

今度は傾いた楕円にノイズがのったようなものが生成されました。

サンプルデータ3: 多重周期性をもつノイズのない系列

もう少し複雑な例として、複数の周期性をもつ系列を考えます。観測値が $y_t = 3 \sin(t) + 2 \cos(t / 3)$ となっている場合を仮定します。可視化したときの見栄えの都合上、1000期分のデータを用意します。

t <- seq(1, 1000, 1)
x1 <- 3 * sin(t) + 2 * cos(t / 3) + rnorm(100, 0, 0.1) 
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

100期までの観測値をプロットすると、以下のようになります。
f:id:tekenuko:20161111001519p:plain
２つの周期的な変動が重なったような振る舞いになっています。

こちらもobsとobs_lagを可視化すると、以下のようになります。
f:id:tekenuko:20161111002017p:plain
単一の周期性をもつ場合に比べると、かなり複雑な構造になりました。線が何回重なるか、何回小さな円っぽい軌道を描くか、などは周期が変わると大幅に変わります。今回は $\sin(t)$ と $\cos(t/3)$ を比較すると、周期が3倍異なるものが２つある、ということで3回回ってもとのところに戻ってくる、という振る舞いになってそうです。

サンプルデータ4: 多重周期性と周期的なスパイクをもつノイズのない系列

さらに複雑な例として、さきほどの系列に周期的なスパイクを加えます。12期ごとに一回、大きな正の値の何かが印加される状況を考えます。

t <- seq(1, 1000, 1)
x1 <- 3 * sin(t) + 2 * cos(t / 3) + ifelse(t %% 12 == 0, 10, 0)
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

100期までの観測値をプロットすると、以下のようになります。
f:id:tekenuko:20161111003001p:plain
たまに大きな値に観測値が振れるような感じになっています。

obsとobs_lagを可視化すると、以下のようになります。
f:id:tekenuko:20161111003023p:plain
何か両端に変なものが湧いてきました。。ちょっとこのままだとわかりづらいので、100000期のデータに増やしてプロットし直してみます。
f:id:tekenuko:20161111003402p:plain
どうやらさきほどの構造のコピーらしきものが両端にできるようです。大きな正の値が印加される、という状況は、ベクトル $\xi_t$ の言葉でいうと、成分 $y_t$ のみが大きな値になる場合と成分 $y_{t-1}$ のみが大きな値になる場合の2パターンが新たに生まれることを意味しています。加えて、印加されるものも周期性を持っているため、新たに生まれた2パターンもさきほどの構造が反映された軌道を描く、といった理解になるんでしょうか。また、こちらの場合、データ点を増やさないと詳細が見えてこなかったということを鑑みると、実データの場合もそうとう蓄積しないとこのような構造が見えない可能性があるな、と思いました。

おわりに

観測値をうまくベクトルにして可視化すると、構造が見える場合があって面白いです。他にも色々なパターンを作れるので、やってみると楽しくて無限に時間を溶かせる感じになります。ここで載せたもの以外にも色々試しており、他にも面白いものもあったりする*1んですが、夜な夜なやりこみすぎると怒られたりするんで、このへんでやめときます。

*1:日中も空いた時間にやって一人でテンション上がってたりしてました（笑）