The parameter estimation problems of the multivariate equation-error systems whose outputs are contaminated by an autoregressive noise process are studied by using the maximum likelihood principle and the coupling identification concept. The key is to decompose the system into several regressive identification subsystems based on the number of the outputs. By means of the negative gradient search, a multivariate coupled maximum likelihood generalized stochastic gradient algorithm is derived. In order to improve the parameter estimation accuracy, a multivariate coupled maximum likelihood multi-innovation generalized stochastic gradient algorithm is proposed by introducing the multi-innovation identification theory. The simulation example indicates that the proposed methods work well.