Global Index (short | long) | Local contents | Local Index (short | long)
[th,it_inf]=bj(z,nn,maxiter,tol,lim,maxsize,Tsamp,trc)
BJ Computes the prediction error estimate of a Box-Jenkins model.
TH = BJ(Z,NN) or TH = BJ(Z,NN,'trace')
TH: Returned as the estimated parameters of the Box-Jenkins model
y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
along with estimated covariances and structure information. For
the exact format of TH, see also THETA.
Z : The output-input data Z=[y u], with y and u being column vectors.
This routine does not work for multi-output systems. Use PEM then.
NN: Initial value and structure information, given either as
NN=[nb nc nd nf nk], the orders and delay of the above model,
or as NN=THI, with THI being a theta matrix of standard format.
In the latter case the criterion minimization is initialized at THI.
For multi-input systems, nb, nf and nk are row vectors, so that nb(k)
nf(k) and nk(k) is the order and delay associated with input # k.
If a *last* argument to oe is entered as 'trace', information about
the criterion minimization is given in the MATLAB Command Window.
Some parameters associated with the algorithm are accessed by
[TH,IT_INF] = BJ(Z,NN,maxiter,tol,lim,maxsize,T,'trace')
See also AUXVAR for an explanation of these and their default values.
See also ARX, ARMAX, IV, N4SID, OE, and PEM.
| This function calls | This function is called by |
|---|---|
function [th,it_inf]=bj(z,nn,maxiter,tol,lim,maxsize,Tsamp,trc)
% L. Ljung 10-1-86,8-27-94
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.5 $ $Date: 1997/12/02 03:42:27 $
if nargin<2
disp('Usage: TH = BJ(Z,NN)')
disp(' TH = BJ(Z,NN,''trace'')')
disp(' TH = BJ(Z,NN,maxiter,tol,lim,maxsize,T,''trace'')')
return
end
% *** Set up default values ***
global XIDiter
gui=0;display=0;stopp=0;flag=[];
eval('flag=findobj(get(0,''children''),''flat'',''tag'',''sitb22'');','');
if ~isempty(flag),if strcmp(get(flag,'vis'),'on')
gui=1;set(XIDiter(7),'userdata',0);
end,end
[Ncap,nz]=size(z);
maxsdef=idmsize(Ncap);
if nargin==8,if strcmp(lower(trc),'trace'),display=1;end,end
if nargin==7,if strcmp(lower(Tsamp),'trace'),display=1;Tsamp=[];end,end
if nargin==6,if strcmp(lower(maxsize),'trace'),display=1;maxsize=[];end,end
if nargin==5,if strcmp(lower(lim),'trace'),display=1;lim=[];end,end
if nargin==4,if strcmp(lower(tol),'trace'),display=1;tol=[];end,end
if nargin==3,if strcmp(lower(maxiter),'trace'),display=1;maxiter=[];end,end
if nargin<7, Tsamp=[];end
if nargin<6, maxsize=[];end
if nargin<5, lim=[];end
if nargin<4, tol=[];, end
if nargin<3, maxiter=[]; end
if display==1,arg='trace';else arg='notrace';end
Tdef=0;if isempty(Tsamp),Tsamp=1;Tdef=1;end
if isempty(maxsize),maxsize=maxsdef;end
if isempty(lim),lim=1.6;end,if isempty(tol),tol=0.01;end
if isempty(maxiter),maxiter=10;end
if Tsamp<0,Tsamp=1;end, if maxsize<0,maxsize=maxsdef;end, if lim<0,lim=1.6;end
if tol<0,tol=0.01;end,if maxiter<0, maxiter=10;end
% *** Do some consistency tests ***
[nr,cc]=size(nn); nu=nz-1;
if nz>Ncap, error('The data should be organized in column vectors')
return,end
if nu==0, error('For time series, use ARMAX instead!'),
return,end
if nr==1
if cc~=3*nu+2,
disp('Incorrect number of orders specified:')
disp('nn should be nn=[nb nc nd nf nk],')
disp('where nb nf and nk are row vectors each of length= # of inputs.')
error(' ')
end
if any(nn<0)
error('All orders must be non-negative.')
end
kz=find(nn(1:nu)==0);
if ~isempty(kz),nn(nu+kz+2)=zeros(1,length(kz));end
nb=nn(1:nu);
if nu>0&sum(nb)==0
nad=nn(nu+2);nc=nn(nu+1);
[th,it_inf]=armax(z(:,1),[nad nc] ,...
maxiter,tol,lim,maxsize,Tsamp,arg);
if isempty(th),return,end
[rth,cth]=size(th);
if cth<6+3*nu,th=[th,zeros(rth,6+3*nu-cth)];end
th(1,3:3+3*nu+3)=[nu zeros(1,nu+1) nn(nu+1:3*nu+2)];
th(3:rth,1:nad+nc)=[th(3:rth,nad+1:nad+nc),th(3:rth,1:nad)];
th(4:3+nad+nc,:)=[th(4+nad:3+nad+nc,:);th(4:3+nad,:)];
th(2,7)=2;
return
end
if nu>1
[th,it_inf]=pem(z,[0 nn(1:3*nu+2)] ,...
[],maxiter,tol,lim,maxsize,Tsamp,arg);
return
end
end % if nr==1
% *** if nn has more than 1 row, it is a theta matrix
% and we jump directly to the minimization ***
%
if nr>1, nu=nn(1,3);
if nu>1,
th=pem(z,nn,[],maxiter,tol,lim,maxsize,Tsamp,arg);
return
end
if nu==0,
error(['This routine is only meaningful for a system with input.',...
' Use ARMAX instead!'])
end
na=nn(1,4);nb=nn(1,5);nc=nn(1,6);nd=nn(1,7);
nf=nn(1,8);nk=nn(1,9);
if na>0,error('na should be zero for this routine!'),return,end
n=nf+nb+nc+nd; t=nn(3,1:n).';ni=max([nd nd+nb+nk-1 nc+nf]);
if nc>0, c=[1 t(nb+1:nb+nc).'];else c=1;,end
if nf>0,f=[1 t(nb+nc+nd+1:n).'];else f=1;,end
if nb>0,b=[zeros(1,nk) t(1:nb).'];,else b=0;,end
if nd>0,d=[1 t(nb+nc+1:nb+nc+nd).'];else d=1;end
w=filter(b,f,z(:,2));v=z(:,1)-w;e=pefilt(d,c,v,zeros(1,ni));
if Tdef,Tsamp=gett(nn);end
end
if nr==1,
%
nb=nn(1); nc=nn(2); nd=nn(3); nf=nn(4);
nk=nn(5); n=sum(nn(1:4));ni=max([nd nd+nb+nk-1 nc+nf]);
%
% *** compute initial estimate via LS and IV ***
%
t=arx(z,[nf nb nk],maxsize,Tsamp,0);
f=fstab([1 t(1:nf).']);
t=iv(z,[nf nb nk],f,[zeros(1,nk) t(nf+1:nf+nb).'],maxsize,Tsamp,0);
if isempty(t)|any(isinf(t))
disp('WARNING: Your input signal is apparently not persistently exciting.')
disp(' Use other input or lower model orders.')
disp(' No model returned.')
th=[];return
end
f=fstab([1 t(1:nf).']);t(1:nf)=f(2:nf+1).';
%
% *** Determine the initial C- and D-estimates by first building a
% high order AR model of the residuals, and then
% using the LS-method ***
%
w=pefilt([zeros(1,nk) t(nf+1:nf+nb).'],f,z(:,2),z(1:ni,1));
v=z(:,1)-w;
c=1; t2=[];
if nc>0
t1=arx(v,n,maxsize,Tsamp,0);
e=pefilt([1 t1.'],1,v);
t2=arx([v e],[nd nc 1],maxsize,Tsamp,0);
%
% ** test the stability of the C-estimate **
%
c=fstab([1 t2(nd+1:nd+nc).']); t2(nd+1:nd+nc)=c(2:nc+1).';
end
%
if nc==0
t2=arx(v,nd,maxsize,Tsamp,0);
end
if nf>0, tf=t(1:nf);else tf=[];end
if nc>0, tc=t2(nd+1:nd+nc);else tc=[];end
if nd>0, td=t2(1:nd);else td=[];end
t=[t(nf+1:nf+nb);tc;td;tf];
%
if nd==0, d=1; else d=[1 t2(1:nd).'];end
e=pefilt(d,c,v,zeros(1,ni));
end
% *** Display initial estimate ***
%
if gui,drawnow,display=get(XIDiter(2),'value');end
Vcap=e'*e/(Ncap-ni);
if display
%clc
disp([' INITIAL ESTIMATE'])
disp(['Current loss:' num2str(Vcap)])
disp(['th-vector'])
disp(t)
end
%
% *** start minimizing ***
%
% ** determine limit for robustification **
if lim~=0, lim=median(abs(e-median(e)))*lim/0.7;end
if lim==0
el=e;
else
[ne,me]=size(e);
ll=ones(ne,me)*lim;la=abs(e)+eps*ll;el=e.*(min(la,ll)./la);
clear ll,clear la
end
testnorm=tol+1; l=0; st=0; nmax=max([nf nb+nk-1 nc nd]);
h=conv(f,c);
%
% ** the minimization loop **
%
while [testnorm>tol l<maxiter st~=1]
l=l+1;
% * compute gradient *
wf=filter(-d,h,w); ef=filter(1,c,e);
vf=filter(-1,c,v); uf=filter(d,h,z(:,2));
M=floor(maxsize/n);
R=zeros(n);Fcap=zeros(n,1);
for k1=nmax:M:Ncap-1
jj=(k1+1:min(Ncap,k1+M));
psi=zeros(length(jj),n);
for k=1:nb, psi(:,k)=uf(jj-k-nk+1);end
for k=1:nc, psi(:,nb+k)=ef(jj-k);end
for k=1:nd, psi(:,k+nb+nc)=vf(jj-k);end
for k=1:nf, psi(:,k+nb+nc+nd)=wf(jj-k);end
R=R+psi'*psi; Fcap=Fcap+psi'*el(jj);
end
if Ncap>M, g=R\Fcap; else g=psi\el(jj);end,grad=Fcap;
testnorm=abs(grad'*g/Vcap/Ncap*100);
%
% * search along the g-direction *
%
[t1,e,el,v,w,V1,c,d,h,st]=...
searchbj(z,t,g,lim,Vcap,nb,nc,nd,nf,nk,ni,display);
if st==1,
[t1,e,el,v,w,V1,c,d,h,st]=searchbj(z,t,grad/trace(R)*length(R),...
lim,Vcap,nb,nc,nd,nf,nk,ni,display);
end
if gui,drawnow,display=get(XIDiter(2),'value');end
if display
%home
disp([' ITERATION # ' int2str(l)])
disp(['Current loss: ' num2str(V1) ' Previous loss: ' num2str(Vcap)])
disp(['Current th prev. th gn-dir'])
disp([t1 t g])
disp(['Norm of gn-vector: ' num2str(norm(g))])
disp(['Expected improvement: ',num2str(testnorm),' %'])
disp(['Achieved improvement: ',num2str((Vcap-V1)*100/Vcap),' %'])
if st==1
disp(['No improvement of the criterion possible along the search direction'])
disp(['Iterations therefore terminated']),end
end %if display
if gui
set(XIDiter(4),'string',int2str(l));
set(XIDiter(5),'string',num2str(V1));
set(XIDiter(9),'string',['(',num2str(Vcap),')']);
set(XIDiter(6),'string',num2str(testnorm));
drawnow
stopp=get(XIDiter(7),'userdata');
end
t=t1; impr=V1-Vcap;Vcap=V1;
if stopp,break,end
end
it_inf=[l,impr,testnorm];
th=zeros(3+n,max([9 n]));
Vcap=e'*e/(length(e)-ni);
if isempty(Vcap)|isempty(t)
disp('WARNING: Your input signal is apparently not persistently exciting.')
disp(' Use other input or lower model orders.')
disp(' No model returned.')
th=[];return
end
th(1,1:9)=[Vcap Tsamp 1 0 nb nc nd nf nk];
th(2,1)=Vcap*(1+n/Ncap)/(1-n/Ncap);
ti=fix(clock); ti(1)=ti(1)/100;
th(2,2:6)=ti(1:5);
th(2,7)=2;
th(3,1:n)=t.';
if maxiter==0,th(2,7)=2;return,end
if Ncap>M, PP=inv(R); else PP=inv(psi'*psi);end
th(4:3+n,1:n)=Vcap*PP;