MATLABで『時系列解析入門』７章

function [aic, Arcoef, s2, maic] = arLevinson(d, M)
%ARLEVINSON
%
% INPUT
% d:data
% M:maximum AR degree (lag)
% 
% OUTPUT
% aic     : AIC
% Arcoef  : AR coefficient (structure)
% s2      : sigma^2
% maic    : The order which provides the minimum AIC

% 前処理
d = d - mean(d);
N = length(d);

% 標本自己共分散を求める
C = xcorr(d, M)/N;
C = C(M+1:end);

% preallocation
aic = zeros(1, M+1);
s2  = zeros(1, M+1);
a   = zeros(M);

% AIC_0
s2(1) = C(1);
aic(1) = N * (log(2 * pi * s2(1))+1) + 2;

% m = 1に対してアルゴリズムを回す
m = 1;
a(m, m) = C(m+1) / s2(1);
s2(2)    = s2(1) * (1 - a(m, m)^2);
aic(m+1) =  N * (log(2 * pi * s2(2)) + 1) + 2*(m + 1);

% m = 2, ... , Mに対してアルゴリズムを回す
for m = 2:M
  % calculate a_m^m
  a(m, m) = (C(m+1) - dot(a(m-1, 1:m-1), C(m:-1:2))) / s2(m);
  % calculate a_i^m
  for i = 1:m-1
     a(m, i) = a(m-1, i) - a(m, m)*a(m-1, m-i); 
  end
  % calculate sigma_{m-1}^2
  s2(m+1) = s2(m) * (1 - a(m, m)^2);
  % calculate AIC_m
  aic(m+1) = N * (log(2 * pi * s2(m+1)) + 1) + 2*(m + 1);
end

% make Arcoef
for m = 1:M
  name = ['order' num2str(m)];
  Arcoef.(name) = a(m, 1:m);
end

% calculate maic
[~, maic] = min(aic);
end

function [aic, Arcoef, s2, maic] = arLeastsquares(d, M)
%ARLEASTSQUARES
%
% INPUT
% d:data
% M:maximum AR degree (lag)
% 
% OUTPUT
% aic     : AIC
% Arcoef  : AR coefficient (structure)
% s2      : sigma^2
% maic    : The order which provides the minimum AIC

% preprocessing
d = d - mean(d);
N = length(d);

% make matrix X Eq.(7.19)
Z = zeros(N-M, M);
for I = 1:M
    Z(:, I) = d(M-I+1:N-I); 
end
y = d(M+1:N);
X = [Z y'];

% Householder(QR decomposition)
[~, R] = qr(X);
S      = R(1:M+1, 1:M+1);

% calculate sigma^2, AIC
s2 = flipud(cumsum(flipud(S(:, end).^2)));
s2 = s2/(N-M);
aic = (N-M) * (log(2*pi*s2)+1) + 2*((0:M)' + 1);

% make Arcoef
for m = 1:M
    name = ['order' num2str(m)];
    res  = S(1:m, 1:m)\S(1:m, M+1);
    Arcoef.(name) = res';
end

% calculate maic
[~, maic] = min(aic);
maic = maic-1;
end

function [aic, V, Arcoef, maic] = mvArLevinson(d, K, M)
%MVARLEVINSON
%
% INPUT
% d:data
% k:dimension of the data
% M:maximum AR degree (lag)
% 
% OUTPUT
% aic     : AIC (size = lag * k)
% Arcoef  : AR coefficient, size = K * K * M * M
% V      : variance-covariance matrix
% maic    : The order which provides the minimum AIC


% preprocessing
d = reshape(d, [], K);
d = d - mean(d);
N = size(d, 1);

% calculate C
C = xcorr(d, M)/N;
C = C(M+1:end, :);
C = reshape(C, M+1, K, K);
C = permute(C,[3 2 1]); 
C0 = C(:, :, 1);
C  = C(:, :, 2:end);

% preallocation
aic = zeros(1, M);
V   = zeros(K, K, M); % = zeros(size(C))
U   = zeros(K, K, M);
A   = zeros(K, K, M, M);
B   = zeros(K, K, M, M);

% Yule Walker
% m = 0
V0   = C0;
U0   = C0;
aic0 = N * (K*log(2*pi) + log(det(V0)) + K) + K*(K+1);

% m = 1 (W_1 = C_1)
A(:, :, 1, 1) = C(:, :, 1) / U0;
B(:, :, 1, 1) = C(:, :, 1)'/ V0;
V(:, :, 1)    = C0 - A(:, :, 1, 1) * C(:, :, 1)';
U(:, :, 1)    = C0 - B(:, :, 1, 1) * C(:, :, 1);
aic(1)        = N * (K*log(2*pi) + log(det(V(:, :, 1))) + K) + K*(K+1) + 2*K^2*1;

% m = 2 ,..., M
for m = 2:M
    W             = C(:, :, m) - sum(pagemtimes(squeeze(A(:, :, m-1, 1:m-1)), C(:, :, m-1:-1:1)), 3);
    A(:, :, m, m) = W /U(:, :, m-1);
    B(:, :, m, m) = W'/V(:, :, m-1);
    for i = 1:m-1
        A(:, :, m, i) = A(:, :, m-1, i) - A(:, :, m, m) * B(:, :, m-1, m-i);
        B(:, :, m, i) = B(:, :, m-1, i) - B(:, :, m, m) * A(:, :, m-1, m-i);
    end
    tC = permute(C, [2 1 3]); 
    V(:, :, m) = C0 - sum(pagemtimes(squeeze(A(:, :, m, 1:m)), tC(:, :, 1:m)), 3);
    U(:, :, m) = C0 - sum(pagemtimes(squeeze(B(:, :, m, 1:m)),  C(:, :, 1:m)), 3);
    aic(m)     = N * (K*log(2*pi) + log(det(V(:, :, m))) + K) + K*(K+1) + 2*K^2*m;
end

% make Arcoef
for m = 1:M
    name = ['order' num2str(m)];
    Arcoef.(name) = squeeze(A(:, :, m, 1:m));
end

% calculate maic
aic = [aic0 aic];
[~, maic] = min(aic);
maic = maic-1;

end