You mean a custom model structure? Read about gray box identification.

Since this is a nonlinear model, you need select an approximate  function to represent the nonlinearity, polynomial function is a good start point. The MA term will cause a big problem, I would suggest use NARX model with many AR terms. Matlab has a built in nlarx function which includes a lot of choices for approximate the nonlinearity.

If you have further question, please let me know.

Cheers

WeiAdd your answer

See my paper below. I don't identify NARMAX model but I use subspace system identification with splines. This might give you a different direction.

In NARMAX systems, since it is in discrete time, all parameters can be properly identified if they enter the system linearly and they only multiplies measured quantities. If linearity does not hold, then it is a research topic that need to be treated case by case.  

If the model structure is known and the model is linear in the parameters, a least squares method can be used to estimate the parameters of the NARX part first and then use the residual as the e(k) data for the estimation of the full NARMAX model. Update the e(k) with the new residuals and repeat the estimation process until a stop criterion is reached.

when the model structure is not known, please refer to the recent book written by Prof. S.A. Billing for the details.

Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains

http://eu.wiley.com/WileyCDA/WileyTitle/productCd-1119943590,subjectCd-EE30.html

If the model structure is already known, as mentioned by Yuzhu, you should use the Extended Least Squares (NARMAX), using the residual e(k) during parameter estimation.

If the model structure is not known, you can use, for instance, Error Reduction Ratio plus some information criteria (Akaike Information Criterion, Bayesian Information Criterion...) to play this role (further details may be found in Billings, S. A. Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains).

If you want to introduce static information during structure selection you may also use the Multiobjective Error Reduction Ratio (see Martins et. al. Improved Structure Detection For Polynomial NARX Models Using a Multiobjective Error Reduction Ratio - http://link.springer.com/article/10.1007%2Fs40313-013-0071-9)

Yes, if the structure of psi is know, you can use the Moore-Penrose pseudoinverse. theta = (psi ' y ) \ ( psi ' psi). If you dont have the parsimonious structure detection method with ERR critérion ready, you could try to fit non-parsimonious structure with regular least square... and increase the orders of the output, input and degree of non-linearities of the polynomial regressors in the psi matrix by "hand". 

There are mutch better approaches tho, im just suggesting a "quick fix" here. I recomend the Nonlinear System identification book by Billings , 2013. Look at chapter 3

what is the structure you used for NARMAX model. Any paper you published in this NARMAC identification

NARMAX toolbox

Matlab code

https://github.com/pkcgritti/matlab-narmaxtool

R code

Code: https://github.com/helonayala/narmax

Paper: https://www.sciencedirect.com/science/article/pii/S2352711019302663