Identification routines want variables to be arranged along the columns, time along rows. Your data "y" is transposed. Also, the order matrix dimensions are wrong (see documentation for the required sizes of orders). I would do something like this:
y = y';
ny = 5; % number of outputs
nu = 1; % number of inputs
data = iddata(y(:,2:6),y(:,1),1);
na = 2*eye(5); % note: na must be ny-by-ny!
nb = 3*ones(5,1); % nb must be ny-by-nu
nc = 0*ones(5,1); % nc must be ny-by-1
nk = zeros(5,1); % nk must be ny-by-nu
model = armax(data,[na nb nc nk])
