時系列データの解析
この例では、timeseries
オブジェクトと関数 regress
を使用して時系列データを可視化し解析する方法を示します。
航空機の乗客データ
まず、1949 年 1 月~1960 年 12 月の期間における航空機の月別乗客数 (千単位) の配列を作成します。
% 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 y = [112 115 145 171 196 204 242 284 315 340 360 417 % Jan 118 126 150 180 196 188 233 277 301 318 342 391 % Feb 132 141 178 193 236 235 267 317 356 362 406 419 % Mar 129 135 163 181 235 227 269 313 348 348 396 461 % Apr 121 125 172 183 229 234 270 318 355 363 420 472 % May 135 149 178 218 243 264 315 374 422 435 472 535 % Jun 148 170 199 230 264 302 364 413 465 491 548 622 % Jul 148 170 199 242 272 293 347 405 467 505 559 606 % Aug 136 158 184 209 237 259 312 355 404 404 463 508 % Sep 119 133 162 191 211 229 274 306 347 359 407 461 % Oct 104 114 146 172 180 203 237 271 305 310 362 390 % Nov 118 140 166 194 201 229 278 306 336 337 405 432 ]; % Dec % Source: % Hyndman, R.J., Time Series Data Library, % http://www-personal.buseco.monash.edu.au/~hyndman/TSDL/. % Copied in October, 2005.
時系列オブジェクトの作成
時系列オブジェクトを作成するときは、データ値と共に時間情報を保持できます。月別データがあるので、日付の配列を作成し、それを Y データと共に使用して時系列オブジェクトを作成します。
yr = repmat((1949:1960),12,1); mo = repmat((1:12)',1,12); time = datestr(datenum(yr(:),mo(:),1)); ts = timeseries(y(:),time,'name','AirlinePassengers'); ts.TimeInfo.Format = 'dd-mmm-yyyy'; tscol = tscollection(ts); plot(ts)
トレンドと季節性の調査
この系列は、線形または 2 次のトレンドがあり、強い季節成分を含んでいるようです。さらに、季節的な変動の大きさは、一般的なレベルの増加と共に増加します。季節的な変動は、対数変換によってさらに一定になる可能性があります。まず、軸のスケールを変更します。
h_gca = gca;
h_gca.YScale = 'log';
季節成分は対数スケールでモデル化した方が簡単のようです。対数変換を使用して新しい時系列を作成します。
tscol = addts(tscol,log(ts.data),'logAirlinePassengers');
logts = tscol.logAirlinePassengers;
次に、月別偏差を重ね合わせて年間平均をプロットします。ここでは、数年の期間にわたる月ごとの変動が一定であるかどうかを判断する必要があります。月-年のフォーマットの行列としてデータを取り扱うこのような操作では、元のデータ行列に対して操作する方が便利です。
t = reshape(datenum(time),12,12); logy = log(y); ymean = repmat(mean(logy),12,1); ydiff = logy - ymean; x = yr + (mo-1)/12; plot(x,ymean,'b-',x,ymean+ydiff,'r-') title('Monthly variation within year') xlabel('Year')
次に、年と月を入れ替えて、年-年のトレンドが各月において一定であるかどうかを判断しましょう。
h_gca = gca; h_gca.Position = [0.13 0.58 0.78 0.34]; subplot(2,1,2); t = reshape(datenum(time),12,12); mmean = repmat(mean(logy,2),1,12); mdiff = logy - mmean; x = mo + (yr-min(yr(:)))/12; plot(x',mmean','b-',x',(mmean+mdiff)','r-') title('Yearly trend within month') xlabel('Month')
トレンドと季節性のモデル化
この系列を、季節成分を含む線形トレンドとしてモデル化してみましょう。
subplot(1,1,1);
X = [dummyvar(mo(:)) logts.time];
[b,bint,resid] = regress(logts.data,X);
tscol = addts(tscol,X*b,'Fit1')
Time Series Collection Object: unnamed Time vector characteristics Start date 01-Jan-1949 End date 01-Dec-1960 Member Time Series Objects: AirlinePassengers logAirlinePassengers Fit1
plot(logts) hold on plot(tscol.Fit1,'Color','r') hold off legend('Data','Fit','location','NW')
このグラフに基づくと、近似は良好のようです。実際のデータと近似値との差は、ここでの目的では問題になるほどの大きさではありません。
ですが、これについてさらに詳しく調べてみましょう。残差は独立しているように見えなければなりません。自己相関 (隣接する残差間の相関) が存在する場合は、それをモデル化し、近似を改善する余地があることが考えられます。残差から時系列を作成し、プロットしてみましょう。
tscol = addts(tscol,resid,'Resid1');
plot(tscol.Resid1)
残差は独立しているようには見えません。むしろ、隣接する残差間の相関は非常に強いように見えます。これは、ダービン-ワトソン検定を使用して正式に検定できます。
[p,dw] = dwtest(tscol.Resid1.data,X)
p = 7.7787e-30
dw = 0.4256
ダービン-ワトソン統計量の低い p 値は、時間と残差が相関していることを示しています。仮説検定の一般的なカットオフは、p<0.05 が有意であるかどうかを判断することです。ここでは非常に小さい p 値が、残差が相関していることを示す強力な証拠となっています。
この自己相関を取り除くために、モデルを変更してみましょう。曲線の一般的な形状は、中央で高く、末端で低くなっています。これは、2 次の項を考慮する必要があることを示しています。ただし、この項を追加しても、自己相関は残るようにも思われます。実際に試してみましょう。
X = [dummyvar(mo(:)) logts.time logts.time.^2];
[b2,bint,resid2] = regress(logts.data,X);
tscol = addts(tscol,resid2,'Resid2');
plot(tscol.Resid2)
[p,dw] = dwtest(tscol.Resid2.data,X)
p = 8.7866e-20
dw = 0.6487
二乗項を追加すると、元の残差プロットで顕著に見られた曲線は取り除かれましたが、プロットと新しいダービン-ワトソン検定の両方において、残差にまだ大きな相関があることが示されています。
このような自己相関は、X 変数ではキャプチャされなかった他の原因によってもたらされている可能性があります。モデルを改善し、相関を削減できる他のデータを収集する必要がありそうです。他のデータが存在しない場合は、単に別のパラメーターをモデルに追加して、自己相関を表すことができます。それを行ってみましょう。ここでは、二乗項を取り除き、誤差に対して自己回帰モデルを使用します。
自己回帰プロセスには次の 2 つの段階があります。
Y(t) = X(t,:)*b + r(t) % regression model for original data r(t) = rho * r(t-1) + u(t) % autoregressive model for residuals
残差系列 r(t)
を一連の個別の値とする場合は、通常の回帰モデルとは異なり、 残差は個別の値で構成される独自の誤差項 u(t)
をもつ自己回帰モデルに従うことができます。
このモデルを作成するために、無名関数 f
を作成して近似値 Yfit
を計算し、Y-Yfit
が u 値をもたらすようにします。
Yfit(t) = rho*Y(t-1) + (X(t,:) - rho*X(t-1,:))*b
この無名関数では、[rho; b]
を単一のパラメーター ベクトル c
に組み合わせます。結果として得られる残差は、相関していない系列にかなり近いものになります。
r = corr(resid(1:end-1),resid(2:end)); % initial guess for rho X = [dummyvar(mo(:)) logts.time]; Y = logts.data; f = @(c,x) [Y(1); c(1)*Y(1:end-1) + (x(2:end,:)- c(1)*x(1:end-1,:))*c(2:end)]; c = nlinfit(X,Y,f,[r; b]); u = Y - f(c,X); tscol = addts(tscol,u,'ResidU'); plot(tscol.ResidU);
まとめ
この例では、Statistics and Machine Learning Toolbox™ の機能を利用して MATLAB® の時系列オブジェクトを使用する方法を説明しています。ts.data
という表記を使用すると、データを抽出し入力として関数に与えることが簡単にできます。また、関数 controlchart
には時系列オブジェクトを直接渡すことができます。
Econometrics Toolbox™ や System Identification Toolbox™ の機能など、時系列用に特別に設計された機能を使用すると、さらに詳細な解析を実行できます。