tfarma_est_mle.m

用于模拟时变非平稳的ARMA过程

function [A, B, exitflag, output]= tfarma_est_mle(y, Aml0, Bml0, re_im, opts)
% function [A, B, exitflag, output]= tfarma_est_mle(y, Aml0, Bml0, re_im, opts)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Exact ML estimator for TFARMA(MAR, LAR; MMA, LMA; N; 1/2)
% models. 
%
%TODO: 
% The switch re_im decides how the parameters are
% optimized! Formulated such that all parameters to be optimized
% are real-valued! (cf. tfar_llfe.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
MAR  =   2;
LAR  =   2;
MMA  =   1;
LMA  =   1;
re_im= 'r';
mo_no= 'm';
tfpm_file_gen;
%-------------
beta = 1/2;

MC= 2*(MAR+MMA);
LC= 2*(LAR+LMA);

Amlref= Aml;
Bmlref= Bml;

y= tfarma_gen(randn(N, 1), Aml, Bml, beta);

qyy= corr_est(y, y, MC, 1/2);
ayy= fft(qyy);
ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
ayy= ayy(N/2+1-3*LC:N/2+1+3*LC, :);

Psi= tf_multiwin(N, MC, LC, 5, 1, 1);
psi= Psi(N/2+1-3*LC:N/2+1+3*LC, N/2+1-MC:N/2+1+MC);

[Aml0, Bml0, INSTAB]= gkmu(ayy.*conj(psi), N, MAR, LAR, MMA, LMA);

opts= optimset('MaxIter', 1000, 'MaxFunEvals', 100000);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions:
[MAR, LAR]= param_dim(Aml0);
[MMA, LMA]= param_dim(Bml0);

% Combine the parameters, do not contain a_{0, l}!
theta_AB= [param_stack_ml(Aml0(:, 2:end)); param_stack_ml(Bml0)];

% Call fminunc:
[X, lX, exitflag, output]= fminunc('tfarma_llfe', theta_AB, opts, y, MAR, LAR, MMA, LMA, re_im);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call fmincon:
[X, lX, exitflag, output]= fmincon('tfarma_llfe', theta_AB, [], [], [], [], -3*ones(size(theta_AB)), 3*ones(size(theta_AB)), [], opts, y, MAR, LAR, MMA, LMA, re_im);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

theta_A= X(1:MAR*(2*LAR+1));
theta_B= X(MAR*(2*LAR+1)+1:end);

Aml= [[zeros(LAR, 1); 1; zeros(LAR, 1)] param_destack_ml(theta_A, MAR-1, LAR)];
Bml= param_destack_ml(theta_B, MMA, LMA);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[param_norm(Aml, Amlref, 0) param_norm(Aml0, Amlref, 0)]
[param_norm(Bml, Bmlref, 1) param_norm(Bml0, Bmlref, 1)]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   2;
LAR  =   2;
MMA  =   0;
LMA  = LAR;
re_im= 'r';
mo_no= 'n';
tfpm_file_gen;
%-------------
beta = 1/2;
B0l= Bml;
MM= 20;
opts= optimset('MaxIter', 300, 'MaxFunEvals', 10000, 'Diagnostics', 'off', 'Display', 'off');

AMLyw= zeros(2*LAR+1, MAR+2, MM);
AMLml= zeros(2*LAR+1, MAR+2, MM);

Iter= [];
NFun= [];
EFLG= [];
for mm= 1:MM
   mm
   y= tfarma_gen(randn(N, 1), Aml, B0l, beta);
   qyy= corr_est(y, y, MAR, 1/2);
   ayy= fft(qyy);
   ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
   ayy= ayy(N/2+1-3*LAR:N/2+1+3*LAR, :);
   [AA, BB]= tfar_est_tfywu(ayy, N);
   Aml0= AA(:, :, end);
   B0l0= BB(:, :, end);
   [A, B, exitflag, output]= tfar_est_mle(y, Aml0, B0l0, re_im, opts);
   Iter= [Iter; output.iterations];
   NFun= [NFun; output.funcCount];
   EFLG= [EFLG; exitflag];
   figure(1);clf;plot(Iter);drawnow
   figure(2);clf;plot(NFun);drawnow
   AMLyw(:, :, mm)= [Aml0 B0l0];
   AMLml(:, :, mm)= [A B];
end;
[mean(Iter) mean(NFun)]
[Myw, Vyw, Byw]= param_mse(AMLyw, [Aml B0l])
[Mml, Vml, Bml]= param_mse(AMLml, [Aml B0l])
[norm(Myw) norm(Mml)]
mean(AMLyw, 3)
mean(AMLml, 3)
sum(EFLG)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%