Main Content

時系列データの解析

この例では、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)

Figure contains an axes object. The axes object with title Time Series Plot:AirlinePassengers, ylabel AirlinePassengers contains an object of type line.

トレンドと季節性の調査

この系列は、線形または 2 次のトレンドがあり、強い季節成分を含んでいるようです。さらに、季節的な変動の大きさは、一般的なレベルの増加と共に増加します。季節的な変動は、対数変換によってさらに一定になる可能性があります。まず、軸のスケールを変更します。

h_gca = gca;
h_gca.YScale = 'log';

Figure contains an axes object. The axes object with title Time Series Plot:AirlinePassengers, ylabel AirlinePassengers contains an object of type line.

季節成分は対数スケールでモデル化した方が簡単のようです。対数変換を使用して新しい時系列を作成します。

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')

Figure contains an axes object. The axes object with title Monthly variation within year, xlabel Year contains 24 objects of type line.

次に、年と月を入れ替えて、年-年のトレンドが各月において一定であるかどうかを判断しましょう。

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')

Figure contains 2 axes objects. Axes object 1 with title Monthly variation within year, xlabel Year contains 24 objects of type line. Axes object 2 with title Yearly trend within month, xlabel Month contains 24 objects of type line.

トレンドと季節性のモデル化

この系列を、季節成分を含む線形トレンドとしてモデル化してみましょう。

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')

Figure contains an axes object. The axes object with title Time Series Plot:logAirlinePassengers, ylabel logAirlinePassengers contains 2 objects of type line. These objects represent Data, Fit.

このグラフに基づくと、近似は良好のようです。実際のデータと近似値との差は、ここでの目的では問題になるほどの大きさではありません。

ですが、これについてさらに詳しく調べてみましょう。残差は独立しているように見えなければなりません。自己相関 (隣接する残差間の相関) が存在する場合は、それをモデル化し、近似を改善する余地があることが考えられます。残差から時系列を作成し、プロットしてみましょう。

tscol = addts(tscol,resid,'Resid1');
plot(tscol.Resid1)

Figure contains an axes object. The axes object with title Time Series Plot:Resid1, ylabel Resid1 contains an object of type line.

残差は独立しているようには見えません。むしろ、隣接する残差間の相関は非常に強いように見えます。これは、ダービン-ワトソン検定を使用して正式に検定できます。

[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)

Figure contains an axes object. The axes object with title Time Series Plot:Resid2, ylabel Resid2 contains an object of type line.

[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);

Figure contains an axes object. The axes object with title Time Series Plot:ResidU, ylabel ResidU contains an object of type line.

まとめ

この例では、Statistics and Machine Learning Toolbox™ の機能を利用して MATLAB® の時系列オブジェクトを使用する方法を説明しています。ts.data という表記を使用すると、データを抽出し入力として関数に与えることが簡単にできます。また、関数 controlchart には時系列オブジェクトを直接渡すことができます。

Econometrics Toolbox™ や System Identification Toolbox™ の機能など、時系列用に特別に設計された機能を使用すると、さらに詳細な解析を実行できます。