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

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

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

% preallocation
s2  = zeros(1, lag);
aic = zeros(1, lag);

% main
for M = 1:lag
    X = zeros(N-M, M);
    for I = 1:M
        X(:, I) = d(M-I+1:N-I);
    end
    y = d(M+1:N);
    theta  = X\y;
    y_est  = X*theta;
    s2(M)  = sum((y - y_est).^2)/(N-M);
    aic(M) = (N-M) * (log(2*pi*s2(M)) + 1) + 2*(M+1);
    
    % save AR coeff
    name = ['order' num2str(M)];
    Arcoef.(name) = theta;
end

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

end

function [aic, V, Arcoef, maic] = mvArNormaleq(d, K, lag)
%MVARNORMALEQ
%
% INPUT
% d   :data
% k   :dimension of the data
% lag :maximum AR degree
% 
% OUTPUT
% aic     : AIC (size = lag * k)
% Arcoef  : AR coefficient
% 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);

% preallocation
V   = zeros(K, K, lag);
aic = zeros(1, lag);

% main
for M = 1:lag
    % make X
    X = zeros(N-M, M*K);
    for I = 1:M
       X(:, K*(I-1)+1:K*I) = d(M-I+1:N-I, :); 
    end
    y = d(M+1:end, :);

    theta = X\y;                % solve normal equation    
    res    = X*theta - y;       % calculate residual
    V(:, :, M)      = cov(res); % save variance covariance matrix
    
    % save AIC
    aic(M) = N * (K*log(2*pi) + log(det(V(:, :, M))) + K) + K*(K+1) + 2*K^2*M;
    
    % save AR coeff
    name = ['order' num2str(M)];
    Arcoef.(name) = permute(reshape(theta, [K, M, K]), [2, 1, 3]);
end

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

function [spec, amp, phase, coh, power, rpower, f] = crossSpectrum(arcoef, v)
% CROSSSPECTRUM 
% 
% INPUTS
% arcoef: AR coeff matrix
% v     : variance-covariance matrix of the noise
% 
% OUTPUTS
% spec  : cross spectrum
% amp   : amplitude spectrum
% phase : phase spectrum
% coh   : coherency
% power : power spectrum decomposition
% rpower: relative power contribution
% f     : frequency

df = 0.001;
f = 0:df:0.5;
N = length(f);
K = size(v, 1);
M = size(arcoef, 1);

% cross spectrum
spec = zeros(K, K, N);

F = zeros(M, N);
for m = 1:M
    F(m, :) = exp(-2*pi*1i*m*f);
end

A = zeros(N, K, K);
for k = 1:K
    for j = 1:K
        for m = 1:M
            A(:, j, k) = A(:, j, k) + (arcoef(m, j, k) * F(m, :))';
        end
        A(:, j, k) = A(:, j, k) -1*(j == k);
    end
end

A = permute(A, [2, 3, 1]);
for I = 1:N
    spec(:, :, I) = A(:, :, I) \ (v / A(:, :, I)');
end
spec = permute(spec, [3, 1, 2]);

% amplitude and phase spectrum
amp = abs(spec);
phase   = angle(spec) + 2*pi;

% coherency
coh = zeros(size(spec));
for k = 1:K
    for j = 1:K
        coh(:, j, k) = amp(:, j, k).^2 ./ (real(spec(:, j, j)) .* real(spec(:, k, k)));
    end
end

% relative power contribution
power = zeros(size(spec));
B = zeros(size(A));
for I = 1:N
    B(:, :, I) = inv(A(:, :, I));
end
B = permute(B, [3, 1, 2]);
for k = 1:K
    for j = 1:K
        power(:, j, k) = v(k, k) * abs(B(:, j, k)).^2;
    end
end

rpower = zeros(size(power));
for k = 1:K
    sumPower = sum(squeeze(power(:, k, :)), 2);
    for j = 1:K
        rpower(:, k, j) = sum(squeeze(power(:, k, 1:j)), 2) ./ sumPower;
    end
end
end