%  arima.m
%  -------------------------------
%     ARIMA 過程の同定
%  -------------------------------
%
N = 240;                        % ある ARIMA過程(ドリフトと年変動を含む) を想定
a = [1 -0.763 0.624];                    %                       A(z)
az = [1 -1]*[a 0;0 a]                    %            (1-z^(-1)) A(z)
azz = [az zeros(1,12-length(az)) -1*az]  % (1-z^(-12))(1-z^(-1)) A(z)
b = [1 0.13 0.02];                       %                       B(z)
%
%                                     % (1-z^(-12))(1-z^(-1)) A(z)*Xn = B(z)*En
%
x = filter(b,azz,5*randn(240,1));     %  サンプル時系列生成    
subplot(421);plot(x); axis([0 240 -1000 1000]);legend('sample Xn',3)
fy = freqz(b,a,N);
subplot(422);semilogy(abs(fy));axis([0 120 1 100]);legend('buried ARMA')
%
y=x(2:N)-x(1:N-1);                     %  階差時系列
subplot(412);plot(y); axis([0 240 -200 200]);legend('difference Yn= Xn+1-Xn',3)
%
z=y(13:N-1)-y(1:N-13);                 %　年変動除去
subplot(413);plot(z); axis([0 240 -100 100]);legend('difference Zn= Yn+12-Yn',3)
%
%  階差時系列を AR過程 で同定
%
ai = aryule(y,20);             % Yule-Walker法  ／次数(20)
fy = freqz(1,ai,N);
subplot(427);semilogy(abs(fy));axis([0 120 1 100]);legend('spectrum of Yn',2)
%
%  年変動除去後の時系列を ARMA過程 で同定
%
ai = aryule(z,30);             % Yule-Walker法  ／次数(30)
hai=impz(1,ai,N);              % インパルス応答
[bs,as] = stmcb(hai,2,2,10)    % 次数(4,4); 10 iterations　（Steiglitz-McBride法）
fz = freqz(bs,as,N);
subplot(428);semilogy(abs(fz));axis([0 120 1 100]);legend('recovered ARMA')
%--------------------------------------------------------
%                                2004.12.20  by K.Tsukada
