MATLABで『時系列解析入門』２章
最終更新：2021/07/14
利用データなど基本的な情報は第一回を見てください.
非専門家が書いています. 十分に批判的に読んで頂くようお願いいたします. 間違い・疑問点などあれば, ぜひコンタクトフォームへ連絡いただけると幸いです.
２章
p.13 時系列の散布図
d = hakusan.velocity;
subplot(221)
plot(d)
xticklabels({}); yticklabels({})
title('方向角速度 y_n')
subplot(222)
histogram(d, 'NumBins', 20)
xticklabels({}); yticklabels({})
title('Histogram of y_n')
subplot(223)
lag = 2;
scatter(d(1:end-lag), d(lag+1:end), 1, '+')
xticklabels({}); yticklabels({})
title('y_{n-2} vs y_n')
subplot(224)
lag = 4;
scatter(d(1:end-lag), d(lag+1:end), 1, '+')
xticklabels({}); yticklabels({})
title('y_{n-4} vs y_n')
元データと, 図2.1 ( p.14 ) の ( a ) , ( c ) , ( e ) に対応している.
p.18 標本自己相関関数
標本平均 $$ \hat{\mu} = \frac{1}{N} \sum_{n = 1}^{N} y_n $$ 標本自己共分散関数 $$ \hat{C_k} = \frac{1}{N} \sum_{n = k+1}^{N} (y_n - \hat{\mu})(y_{n-k} - \hat{\mu}) % = \mathbb{E}[(y_n - \hat{\mu})(y_{n-k} - \hat{\mu})] $$ 標本自己相関関数 $$ \hat{R_k} = \frac{\hat{C_k}}{\hat{C_0}} $$
今回求めるのは, 標本自己相関関数である.
MATLABには, xcorrという関数があるが, これが求めてくれるのは, documentを見る限り,
$$
\sum_{n = k+1}^{N} y_n y_{n-k}
% \mathbb{E}[y_n y_{n-k}]
$$
なので, 先に $\hat{\mu}$ を引いておく必要がある.
さらに, $\hat{C_k}$ を求めるためには, $N$ で割る必要もあるだろう.
また, 標本自己相関関数を求めるためには, 標本自己共分散関数を $\hat{C_0}$ でわってもいいし, optionで'normalized'をつけてもいい.
subplot(211)
d = hakusan.velocity;
d = d - mean(d);
[c, lags] = xcorr(d, 50, 'normalized');
bar(lags(51:end), c(51:end));
ylim([-1.0 1.0])
ylabel('\it{R_k}')
subplot(212)
d = whard;
d = d - mean(d);
[c, lags] = xcorr(d, 50, 'normalized');
bar(lags(51:end), c(51:end));
ylim([-1.0 1.0])
ylabel('\it{R_k}')
自分のデータにもつかってみた.
12データ周期くらいのサイクルがある. 月単位の気候データだとばれそうである.
p.23 標本相互相関関数
今度は多変量. 前節の式に, どの時系列かを示す $i$, $j$ が入るだけである.
標本平均 ( ベクトル ) $$ \hat{\mu}(i) = \frac{1}{N} \sum_{n = 1}^{N} y_n(i) $$ 標本相互共分散関数 $$ \hat{C_k}(i, j) = \frac{1}{N} \sum_{n = k+1}^{N} (y_n(i) - \hat{\mu}(i))(y_{n-k}(j) - \hat{\mu}(j)) $$ 標本相互共分散関数 $$ \hat{R_k}(i, j) = \frac{\hat{C_k}(i, j)}{\sqrt{{\hat{C_0}(i, i)}{\hat{C_0}(j, j)}}} $$
標本相互共分散関数, 標本相互共分散関数は $(i, j)$ がついて複雑だが, 自己のときは「自分 ( の系列 ) と自分」の組み合わせしかなったのが, 「自分と他人」の組み合わせに増えただけである.
MATLABの xcorr はもともと, 相互相関を求めるための関数なので, 簡単に実践できる.
yoko = hakusan.yoko;
tate = hakusan.tate;
dakaku = hakusan.dakaku;
yoko = yoko - mean(yoko);
tate = tate - mean(tate);
dakaku = dakaku - mean(dakaku);
% 対角成分, つまり標本自己相関関数
c = xcorr(yoko, 50, 'normalized');
subplot(331); plot(c(51:end));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({});
c = xcorr(tate, 50, 'normalized');
subplot(335); plot(c(51:end));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
c = xcorr(dakaku, 50, 'normalized');
subplot(339); plot(c(51:end));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
% 非対角成分, つまり標本相互相関関数
c = xcorr(yoko, tate, 50, 'normalized'); % 横揺れ vs 縦揺れ
subplot(332); plot(c(51:end));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
subplot(334); plot(c(51:-1:1));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
c = xcorr(yoko, dakaku, 50, 'normalized'); % 横揺れ vs 舵角
subplot(333); plot(c(51:end));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
subplot(337); plot(c(51:-1:1));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
c = xcorr(tate, dakaku, 50, 'normalized'); % 縦揺れ vs 舵角
subplot(336); plot( c(51:end));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
subplot(338); plot(c(51:-1:1));
ylim([-1.0 1.0]); xlim([0 50]); yline(0.0); xticklabels({}); yticklabels({})
主な参考文献
