The ARRMarkov-least squares method is extended to multivariable systems. This method explicitly determines the Markov parameters (impulse response coefficients) of a process using process input-output data and a standard least-squares (LS) algorithm. The parameter estimates are consistent and have tighter confidence bounds than those produced by other linear regression methods. The Interactor matrix, which defines the time delays in multivariable systems, can be directly estimated from the Markov parameters. From simulation results it is observed that the Markov parameters estimated by the ARMarkov-LS method an the closest to the actual Markov parameters irrespective of the system order and lead to a better estimate of the interactor matrix than other linear regression methods such as Correlation Analysis, ARX, FIR, etc. The identified Markov parameters and/or the time-delay/interactor matrix can be used directly in the design of model predictive controllers and control loop performance assessment.