What's wrong with just doing sys = armax(data,[na,nb,nc,nk]);
as per the Mathworks manual? sys is the variable. You could do sys[i] = ... to store several of them in an array (or perhaps sys{i}).
Then, to compare them, you could use getcov()
and compare the covariances of each model. These covariances tell you how "uncertain" the model is.
Another approach might be to look at the model residuals. If the residuals are white (uncorrelated), then the model is good. The relatively energy / power of the residuals will be a measure of how accurate the model matches the data (the smaller the residual, the better the model).