function [XX] = MakeRegressionMatrix(X0,Memory,order)
% Making input to a direct order identification of Volterra
%
% Input = X0(1xN)
% Memory = memory term of the series, Integer
% Order = Order of the model, Integer
%Outputs:
%XX = Full regression matrix structure
 
 
% Make first order kernel input
X1 = zeros(length(X0),Memory+1);

if isrow(X0)
    X0 = X0';
end


%% Make first order Kernels
X1(:,1) = X0;

for i=1:1:Memory
X1(1+i:end,i+1)=X0(1:end-i);
end

XX = X1;
X2 = NaN;
X3 = NaN;


%% Make 2nd order kernels
if order > 1
    
    X2a = zeros(length(X0),(Memory+1)^2);
    
    X2b = repmat(X1,1,Memory+1); 
    X2aa = repmat(X0,[1 (Memory+1)]);
    
    % 2nd order direct kernel for X
    for i=1:1:(Memory+1)
       X2b(:,(i-1)*Memory+i:i*Memory+i)=X1;
       X2a(i:end,(i-1)*Memory+i:i*Memory+i)=X2aa(1:end-i+1,:);
    end
     
    X2 = X2a.*X2b;
    XX = [XX X2];
end
    
    %% Make 3rd order kernels
if order > 2
    
    X3aa = repmat(X2,1,(Memory+1));

    X3bb = zeros(size(X3aa));
    X3b1 = repmat(X2aa,1,Memory+1); 

    for i=1:1:(Memory+1)
       X3bb(i:end,(i-1)*(Memory+1)^2+1:i*(Memory+1)^2)=X3b1(1:end-i+1,:);
    end

    X3 = X3bb.*X3aa;
    XX = [XX X3];
end

if order > 3
    error('4th-order Volterra og higher is not implemented')

end

end

