Before fitting a model type (ARMAX, OE...), I normally recommend using matlab tools to estimate the delay(s) between your input(s) and output, before identifying dynamic parameters. I am assuming your model is linear and time-invariant) Options are:

1)

Use time domain approaches like cross-correlation to get an estimate of system delays if your data record is continuous (time invariant) and sufficiently long. This may be a noisy estimate difficult to locate on a noisy correlation signal, but it should at least define a reasonable range for the delay. AND/OR

Use a 4th order ARX formulation (fixed complexity high enough to overmodel your system. This can work well even on shorter data sets found in hybrid or time-varying systems. Here you fix the range and resolution steps for the delay search;  estimate the parameters for each potential delay. Pick the delay that minimizes the MSE on a validation data set. The parameters here are not preserved, because you are only trying to refine the delay value using an over-complex model.

2) Now fix the delay  to the estimated value, and find the approximate model order required with the ARX form by comparing QF (BIC, MDL, AIC...) and MSE on a validation data set. Pick the best ARX form using in addition other factors like whiteness and zero-mean residuals.

3) In the previous step  the ARX order may still be too large for your system, so plot the pole-zeroes. Matlab allows plotting estimated pole-zeroes with their confidence intervals, which allows you to detect any pole-zero cancellations (confidence ellipse overlap). Removing these provides a trimmed model order at least for ARX form.

5) The surviving number of poles and zeroes will determine the required maximum complexity for polynomials of any other model form (OE, ARMAX. BJ) for the input processes, while first or second order should suffice for the noise processes.

Now you attempt true identification, with the delay fixed to a value found above, fixed model complexity (numbers poles and zeroes), and search for the optimal complexity for each model form, using usual criteria of white residuals with zero mean, model description length (MDL), MSE etc . Parsimony normally wins: the models with the same criteria levels, but fewer parameters are preferred. However you may have different goals, such as in validation data set: matching responses to special inputs, intended use of the model which may highlight some parameters, or a desire for especially robust estimates of key parameters.

I recommend that any tools you develop like this for identification be tested on simulated data of similar characteristics (but known parameters = gold standard) to develop a feeling for the performance and noise tolerance of your problem.

You can try different criteria such as the AIC and BIC. Other options include genetic algorithm and the OFR algorithm (Orthogonal Forward Regression algorithms).

A cross-correlation based model validation can also solve your problem.

Thank you for your helpful answer. I once thought using genetic algorithm, but It's too slow; therefore, I will try your other suggestions. 

I strongly suggest using Young Information Criterion (YIC) instead of AIC or BIC ones. it includes the covariance matrix and give you more reliable results.

Before fitting a model type (ARMAX, OE...), I normally recommend using matlab tools to estimate the delay(s) between your input(s) and output, before identifying dynamic parameters. I am assuming your model is linear and time-invariant) Options are:

1)

Use time domain approaches like cross-correlation to get an estimate of system delays if your data record is continuous (time invariant) and sufficiently long. This may be a noisy estimate difficult to locate on a noisy correlation signal, but it should at least define a reasonable range for the delay. AND/OR

Use a 4th order ARX formulation (fixed complexity high enough to overmodel your system. This can work well even on shorter data sets found in hybrid or time-varying systems. Here you fix the range and resolution steps for the delay search;  estimate the parameters for each potential delay. Pick the delay that minimizes the MSE on a validation data set. The parameters here are not preserved, because you are only trying to refine the delay value using an over-complex model.

2) Now fix the delay  to the estimated value, and find the approximate model order required with the ARX form by comparing QF (BIC, MDL, AIC...) and MSE on a validation data set. Pick the best ARX form using in addition other factors like whiteness and zero-mean residuals.

3) In the previous step  the ARX order may still be too large for your system, so plot the pole-zeroes. Matlab allows plotting estimated pole-zeroes with their confidence intervals, which allows you to detect any pole-zero cancellations (confidence ellipse overlap). Removing these provides a trimmed model order at least for ARX form.

5) The surviving number of poles and zeroes will determine the required maximum complexity for polynomials of any other model form (OE, ARMAX. BJ) for the input processes, while first or second order should suffice for the noise processes.

Now you attempt true identification, with the delay fixed to a value found above, fixed model complexity (numbers poles and zeroes), and search for the optimal complexity for each model form, using usual criteria of white residuals with zero mean, model description length (MDL), MSE etc . Parsimony normally wins: the models with the same criteria levels, but fewer parameters are preferred. However you may have different goals, such as in validation data set: matching responses to special inputs, intended use of the model which may highlight some parameters, or a desire for especially robust estimates of key parameters.

I recommend that any tools you develop like this for identification be tested on simulated data of similar characteristics (but known parameters = gold standard) to develop a feeling for the performance and noise tolerance of your problem.

Thank you Professor for your very detailed response which will change my approach to linear SI. Before I wrote this question, I tried the Matlab's "Infer Input Delay" function under Nonlinear Identification tab; however, the delays that I found from there were inconsistent. Actually, my system is nonlinear, but I want to see the potential of linear identification methods on my system. 

I suggest that for a family of model candidates the following criterions are tested: Information criterions (AIC, BIC), Error performance measures (MAE, MAPE, RMSE,...), Error characteristics (no significant autocorrelation, Error is approximately white noise, JB test, etc), Homoskedastisity tests, Multicollinearity tests for input predictors, t tests for parameters significance, overall F test, and so on. Then you choose the most suitable model (lags, structure,...), which fullfils the majority of mentioned criterions.