We consider a parameter estimation problem in a general class of stochastic multiple-inputs multiple-outputs Wiener models, where the likelihood function is, in general, analytically intractable. When the output signal is a scalar independent stochastic process, the likelihood function of the parameters is given by a product of scalar integrals. In this case, numerical integration may be efficiently used to approximately solve the maximum likelihood problem. Otherwise, the likelihood function is given by a challenging multidimensional integral. In this contribution, we argue that by ignoring the temporal and spatial dependence of the stochastic disturbances, a computationally attractive estimator based on a suboptimal predictor can be constructed by evaluating scalar integrals regardless of the number of outputs. Under some conditions, the convergence of the resulting estimators can be established and consistency is achieved under certain identifiability hypothesis. We highlight the relationship between the resulting estimators and a recently proposed prediction error method estimator. We also remark that the method can be used for a wider class of stochastic nonlinear models. The performance of the method is demonstrated by a numerical simulation example using a 2-inputs 2-outputs model with 9 parameters.