We present a re-parameterization of vector autoregressive moving average (VARMA) models that allows estimation of parameters under the constraints of causality and invertibility. The parameter constraints associated with a causal invertible VARMA model are highly complex. Currently there are no procedures that can maintain the constraints in the estimated VARMA process, except in the special case of a vector autoregression (VAR), where some moment based causal estimators are available. Even in the VAR case, the available likelihood based estimators are not causal. The maximum likelihood estimator based on the full likelihood that does not condition on the initial observations by definition satisfies the causal invertible constraints but optimization of the likelihood under the complex constraints is an intractable problem. The commonly used Bayesian procedure for VAR often has posterior mass outside the causal set because the priors are not constrained to the causal set of parameters. We provide an exact mathematical solution to this problem. An $m$-variate VARMA$(p, q)$ process contains $(p+ q) m^2 + \binom{m+1}{2}$ parameters, which must be constrained to a subset of Euclidean space in order to guarantee causality and invertibility. This space is implicitly described in this paper, through the device of parameterizing the entire space of block Toeplitz matrices in terms of positive definite matrices and orthogonal matrices. The parameterization has connection to Schur- stability of polynomials and the associated Stein transformation that are often used in dynamical systems literature. As an important by-product of our investigation, we generalize a classical result in dynamical systems to provide a characterization of Schur stable matrix polynomials.