tfarma_est_mlu1.m

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

function [A, B, exitflag, output]= tfarma_est_mlu1(Ayy, Aml0, Bml0, re_im, opts)
% function [A, B, exitflag, output]= tfarma_est_mlu1(Ayy, Aml0, Bml0, re_im, opts)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Asymptotic ML estimator for TFARMA(MAR, LAR; MMA, LMA; N; 1/2)
% models. The switch re_im decides how the parameters are
% optimized! Crucial? 

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

%Bml(LMA+1, 1)= 2;
Amlref= Aml;
Bmlref= Bml;

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

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

opts= optimset('MaxIter', 50, 'MaxFunEvals', 100000, 'GradObj', 'on', 'DerivativeCheck', 'off', 'Display', 'iter', 'Diagnostics', 'on');

Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
Ayy= Aalpha;

[Aml0, Bml0, INSTAB]= gkmu(Ayy(N/2-3*LC+1:N/2+3*LC+1, N/2-MC+1:N/2+MC+1), N, MAR, LAR, MMA, LMA);
param_norm(Aml0, Aml, 1)
param_norm(Bml0, Bml, 1)


y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
%plot(real(y))
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, :).*conj(psi);

[Aml0, Bml0, INSTAB]= gkmu(ayy, N, MAR, LAR, MMA, LMA);
param_norm(Aml0, Aml, 1)
param_norm(Bml0, Bml, 1)

%Aml0= Aml;
%Aml0(:, 2:end)= Aml0(:, 2:end)+.01*(randn(size(Aml(:, 2:end))) + j*randn(size(Aml(:, 2:end))));
%Bml0= Bml;
%Bml0(:, 2:end)= Bml0(:, 2:end)+.01*(randn(size(Bml(:, 2:end))) + j*randn(size(Bml(:, 2:end))));

qyy= corr_est(y, y, -1, 1/2);
ayy= fft(qyy);
ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
Ayy= ayy.*conj(Psi);

%Ayy= tfarma_ambi(Aml, Bml, N, 1/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)];

theta_ABreim= [real(param_stack_ml(Aml0(:, 2:end))); ...
	       imag(param_stack_ml(Aml0(:, 2:end))); ...
	       real(param_stack_ml(Bml0)); ...
	       imag(param_stack_ml(Bml0))];

Pyy= real(nm_to_nk(ml_to_nm(Ayy)));

% Call fminunc:
%[X, lX, exitflag, output]= fminunc('tfarma_llfu', theta_AB, opts, Pyy, MAR, LAR, MMA, LMA, re_im);
[X, lX, exitflag, output]= fminunc('tfarma_llfu1', theta_ABreim, opts, Pyy, MAR, LAR, MMA, LMA, re_im);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tfarma_llfu(theta_AB, Pyy, MAR, LAR, MMA, LMA, re_im);

AB= reshape(theta_AB, 2*LAR+1, length(theta_AB)/(2*LAR+1));
tfar_llfu(AB, Pyy, re_im)

% Call fmincon:
[X, lX, exitflag, output]= fmincon('tfarma_llfu', theta_AB, [], [], [], [], -3*ones(size(theta_AB)), 3*ones(size(theta_AB)), [], opts, Pyy, MAR, LAR, MMA, LMA, re_im);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X= [X(1:MAR*(2*LAR+1)) + j*X(MAR*(2* LAR+1)+1:2*MAR*(2*LAR+1));
    X(2*MAR*(2*LAR+1)+1:2*MAR*(2*LAR+1)+(MMA+1)*(2*LMA+1)) ...
    + j*X(2*MAR*(2*LAR+1)+(MMA+1)*(2*LMA+1)+1:end)];

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[param_norm(A, Aml, 1) param_norm(Aml0, Aml, 1)]
[param_norm(B, Bml, 1) param_norm(Bml0, Bml, 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    =  64;
MAR  =   2;
LAR  =   2;
MMA  =   1;
LMA  =   1;
re_im= 'r';
mo_no= 'm';   
tfpm_file_gen;
%-------------
beta = 1/2;
MM= 50;

AMLgk= zeros(2*LAR+1, MAR+1, MM);
AMLml= zeros(2*LAR+1, MAR+1, MM);
BMLgk= zeros(2*LMA+1, MMA+1, MM);
BMLml= zeros(2*LMA+1, MMA+1, MM);

Iter= [];
NFun= [];
EFLG= [];


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

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);
psi= ones(size(psi));

opts= optimset('MaxIter', 50, 'GradObj', 'on', 'DerivativeCheck', 'off', 'Display', 'iter', 'Diagnostics', 'on');

Psi= tf_multiwin(N, MAR, LAR, 5, 1, 1);
for mm= 1:MM
   mm

   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, :).*conj(psi);
   Qyy= corr_est(y, y, -1, 1/2);
   Ayy= fft(Qyy);
   Ayy= [Ayy(N/2+1:N, :); Ayy(1:N/2, :)].*conj(Psi);

%   [Aml0, Bml0, INSTAB]= gkmu(ayy, N, MAR, LAR, MMA, LMA);
%   param_norm(Aml0, Aml, 1)
%   param_norm(Bml0, Bml, 1)

   Aml0= Aml;
   Bml0= Bml;
   Aml0(:, 2:end)= Aml0(:, 2:end) + (randn(size(Aml0(:, 2:end))) + j*randn(size(Aml0(:, 2:end))))/20;
   Bml0= Bml0  + (randn(size(Bml0)) + j*randn(size(Bml0)))/20;
   
   [A, B, exitflag, output]= tfarma_est_mlu1(Ayy, Aml0, Bml0, 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
   AMLgk(:, :, mm)= Aml0;
   BMLgk(:, :, mm)= Bml0;
   AMLml(:, :, mm)= A;
   BMLml(:, :, mm)= B;
end;
[mean(Iter) mean(NFun)]
[MAgk, VAgk, BAgk]= param_mse(AMLgk, Aml);
[MBgk, VBgk, BBgk]= param_mse(BMLgk, Bml);
[MAml, VAml, BAml]= param_mse(AMLml, Aml);
[MBml, VBml, BBml]= param_mse(BMLml, Bml);

[norm(MAgk) norm(MAml)]
[norm(MBgk) norm(MBml)]
[norm(VAgk) norm(VAml)]
[norm(VBgk) norm(VBml)]
[norm(BAgk) norm(BAml)]
[norm(BBgk) norm(BBml)]
%mean(AMLgk, 3)
mean(AMLml, 3)
%mean(BMLgk, 3)
mean(BMLml, 3)
sum(EFLG)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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    =  64;
MAR  =   2;
LAR  =   2;
MMA  =   1;
LMA  =   1;
re_im= 'r';
mo_no= 'm';   
tfpm_file_gen;
%-------------
beta = 1/2;
alpha= 1/2;

theta_AB= [param_stack_ml(Aml(:, 2:end)); param_stack_ml(Bml)];
theta_ABreim= [real(param_stack_ml(Aml(:, 2:end))); ...
	       imag(param_stack_ml(Aml(:, 2:end))); ...
	       real(param_stack_ml(Bml)); ...
	       imag(param_stack_ml(Bml))];

Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
Ayy= Aalpha;
Pyy= real(nm_to_nk(ml_to_nm(Ayy)));




[l, v]= tfarma_llfu1(theta_ABreim, Pyy, MAR, LAR, MMA, LMA, re_im);
delt= zeros(size(theta_ABreim));
delt(1)= .001* theta_ABreim(1);
[l1, v1]= tfarma_llfu1(theta_ABreim+delt, Pyy, MAR, LAR, MMA, LMA, re_im);
(l1-l)/delt(1)


ff= 4

[l, v]= tfarma_llfu(theta_AB, Pyy, MAR, LAR, MMA, LMA, re_im);
delt= zeros(size(theta_AB));
delt(ff)= .04* theta_AB(ff);
[l1, v1]= tfarma_llfu(theta_AB+delt, Pyy, MAR, LAR, MMA, LMA, re_im);
(l1-l)/delt(ff)
-1/2*v1(ff)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%