function [mbest,minaic,pbest,qbest,Pbest,Qbest]=sarmabat(x,pvec,qvec,Pvec,Qvec,T)
%
%   [mbest,minaic,pbest,qbest,Pbest,Qbest]=sarmabat(x,pvec,qvec,Pvec,Qvec,T)
%
% This functions computes SARMA or multiplicative (p,q) x (P,Q) models for 
% (p,q,P,Q) in (pvec x qvec x Pvec x Qvec); it returns the best according to AIC
% where AIC has been modified to account for fixed parameters
%
% x = input data
% pvec = vector of p's ; set pvec=[0] for no AR
% qvec = vector of q's; set qvec=0[] for no MA
% Pvec = vector of P's
% Qvec = vector of Q's
% T period for multiplicative model

% now estimate ARMA model; ARMAX is AX=BU + Ce so identify phi with A and theta with C
nx=length(x);
np=length(pvec);
nP=length(Pvec);
nq=length(qvec);
nQ=length(Qvec);
aicsave=-99*ones(np,nP,nq,nQ);
fpesave=-99*ones(np,nP,nq,nQ);
minaic=1e+6;
for pp=1:np     
    p=pvec(pp);
    for PP=1:nP
        P=Pvec(PP);
        for qq=1:nq
            q=qvec(qq);
            for QQ=1:nQ
                Q=Qvec(QQ);
                if p+P+q+Q ~=0
                    %  now we must set up an initial model structure and freeze the correct coeffs
                    % to do this we will first use dummy coefficients of 1 everywhere
                    phi=ones(1,p+1);      % phi becomes A for matlab 
                    Phi=ones(1,P+1);
                    PhiT=zeros(1,1+(length(Phi)-1)*T);
                    indexes=[1:T:length(PhiT)];
                    PhiT(indexes)=Phi;  % PhiT now has sparse ones
                    
                    theta=ones(1,q+1);    % theta becomes C for matlab
                    Theta=ones(1,Q+1);
                    ThetaT=zeros(1,1+(length(Theta)-1)*T);
                    indexes=[1:T:length(ThetaT)];
                    ThetaT(indexes)=Theta;
                    
                    % now build the product polynomials
                    thetaM=conv(theta,ThetaT);
                    phiM=conv(phi,PhiT);
                    
                    % this just changes to matlab notation to help programming
                    A=phiM; % starting poly
                    Bdum=[];
                    C=thetaM; % starting poly
                    Ddum=[];
                    
                    mi=idpoly(A,Bdum,C,Ddum); % A is phi and C is theta here
                    
                    zphiM=find(mi.a==0);    % addresses where phiM==0
                    pset=[];                % pset is the collection of parameter numbers to be frozen
                    if length(zphiM)     
                        zphiM=zphiM-1;      % adjust down as leading 1 does not count
                        pset=[pset zphiM];  
                    end
                    zthetaM=find(mi.c==0);  % addresses where thetaM==0
                    if length(zthetaM)
                        zthetaM=zthetaM-1+mi.na;    % subtract 1 as leading 1 does not count and
                        pset=[pset zthetaM];        % add the number of parameters from A=phiM    
                    end                             % that can be changed (the a part of mi.par=length(mi.a)-1)
                    % now set the starting parameters to zero except the leading ones
                    A=zeros(size(phiM));
                    A(1)=1;
                    C=zeros(size(thetaM));
                    C(1)=1;
                    mi=idpoly(A,Bdum,C,Ddum);  % set mi again 
                    set(mi,'FixedParameter',pset);   % set the fixed places by number
                    
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% now finally ready to do armax
                    
                    m=armax(x,mi);          % m is a structure with the model stuff in it
                    resids=pe(m,x);         % pe returns the prediction errors
                    nres=length(resids);
                    rhores=acf(resids,1); % this returns the acf normalized by the variance
                    nrho=length(rhores);             
                    % next compute the Ljung - Box statistic and P-values
                    deltak=floor(nrho/10)+1;
                    kvec=[p+q+P*T+Q*T+1:deltak:p+q+P*T+Q*T+1+4*deltak];
                    for kk=1:5
                        Qsum=0;
                        for j=2:kvec(kk)+1
                            Qsum=Qsum+(rhores(j).^2)/(nx-j);
                        end
                        Qsum=nx*(nx-1)*Qsum;
                        ljpv(kk)=1-chi2cdf(Qsum,kvec(kk)-p-q-P-Q);  % df=kvec(kk)-less no. pars
                    end
                    
                    fpeval=fpe(m);
                    fpesave(pp,PP,qq,QQ)=fpeval;

                    nval = m.EstimationInfo.DataLength;
                    vval = m.EstimationInfo.LossFcn;
                    totalpars=m.na+m.nc;  % BD adds 1 here, probably for the mean; we omit for now
                    netpars=totalpars-length(m.fix);
                    aicval=log(vval)+2*netpars/nval;
                    % aicval=aic(m); % the standard way not taking fixed ones into account
                    aicsave(pp,PP,qq,QQ)=aicval;
                    bicval=log(vval)+netpars*log(nval)/nval;
                    aiccval=log(vval)+2*netpars/(nval-2*netpars); % my best understanding of adapting
                    if aicval < minaic                            % BD to Choi and matlab
                        minaic=aicval; % save the min
                        pbest=p;
                        qbest=q;
                        Pbest=P;
                        Qbest=Q;
                        mbest=m;
                        % put this here to only prints the ones that improve
                        %disp(sprintf('(p,q)(P,Q)=(%d,%d)(%d,%d) aic=%g  fpe=%g',p,q,P,Q,aicval,fpeval));
                    end
                    % next is routine print
                    disp(sprintf('(p,q)(P,Q)=(%d,%d)(%d,%d) aic=%g  aicc=%g bic=%g fpe=%g',...
                        p,q,P,Q,aicval,aiccval,bicval,fpeval));

                    % next was for debugging
                    %disp(sprintf('(p,d,q)(P,D,Q)=(%d,%d,%d)(%d,%d,%d) aic=%g  fpe=%g',p,d,q,P,D,Q,...
                     %   aicsave(pp,PP,qq,QQ),fpesave(pp,PP,qq,QQ)));
                    LJPV=[kvec;ljpv];
                    % next is routine print
                    %disp(sprintf('Ljung-Box P-values : '));
                    %disp(sprintf('  K=%d P-v=%6.4f \n',LJPV(:)));
                    

                end % if p+P+q+Q 
            end % QQ loop
        end     % qq loop
    end         % PP loop
end             % pp loop