

%y outputs (n+p)x1
%U inputs (n+p)x m, m=number of inputs
%p = number of predictor coefficients to estimate
%lab = "yy", "ny" or "yn" for, respectively, ARMAX, Output error, Time series
%n = (dimension of y) - p, data in y^+

function [A,m,nv,yred,ytrue]=Abuild(y,U,p,lab);

fat=2/3;%establishes ARX order
ARXorder=p;%maximum ARX order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% computing A,T,m,yred, ytrue %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(y,1)-p;
m1=size(y,2);%number of outputs
m2=size(U,2);%number of inputs
T=[];
if lab=='yy'%ARMAX
    m=m1+m2;
    for i=1:m1;
        for j=1:n;
            f=y(:,i)';
            A(j,:,i)=fliplr(f([j:(j+p-1)]));
        end
        T=[T A(:,:,i)];
    end
    for i=(m1+1):m;
        for j=1:n;
            f=U(:,i-m1)';
            A(j,:,i)=fliplr(f([j:(j+p-1)]));
        end
        T=[T A(:,:,i)];
    end
elseif lab=='ny';%output error
    m=m2;
    for i=1:m2;
        for j=1:n;
            f=U(:,i)';
            A(j,:,i)=fliplr(f([j:(j+p-1)]));
        end
        T=[T A(:,:,i)];
    end
elseif lab=='yn';%time series
    m=m1;
    for i=1:m1;
        for j=1:n;
            f=y(:,i)';
            A(j,:,i)=fliplr(f([j:(j+p-1)]));
        end
        T=[T A(:,:,i)];
    end
end
yred=y(p+1:n+p,:);
ytrue=yred;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% computing nv (innovation variance) %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pr=min(ARXorder,floor((n*fat)/m));
pr=min(p,pr);
Tr=[];
for i=1:m;
    Tr=[Tr A(:,1:pr,i)];
end
for i=1:m1;
    [X] = lscov(Tr,yred(:,i));
    SD=sqrt(sum((yred(:,i)-Tr*X).^2)/(n-size(Tr,2)));
    nv(i)=SD^2;
end
%q=pem([y U]);nv=q.NoiseVariance;
nv=mean(nv);
SD=sqrt(nv);

  

