HELP: how to estimate coefficient matrices of VARMA model

Hi eveyone, I am a freash man to use SAS. In my work, there are a batch of observations of (demand, price), because of their correlation, I want to use vector ARMA model to predict the future demand and price. But I do not know how do use SAS to estimate the coefficient matrices and the covariance matrix of white noise of the general VARMA model
as below

(I-Phi(1)B-Phi(2)B^2-...-Phi(p)B^p)Z(t)=(I-Theta(1)B-Theta(2)B^2-...-Theta(p)B^q)a(t)
where Z(t) and a(t) are vectors, a(t) is the noise?