  We construct the general permutation invariant Gaussian 2-matrix model for
matrices of arbitrary size $D$. The parameters of the model are given in terms
of variables defined using the representation theory of the symmetric group
$S_D$. A correspondence is established between the permutation invariant
polynomial functions of the matrix variables (the observables of the model) and
directed colored graphs, which sheds light on stability properties in the large
$D$ counting of these invariants. The refined counting of the graphs is given
in terms of double cosets involving permutation groups defined by the local
structure of the graphs. Linear and quadratic observables are transformed to an
$S_D$ representation theoretic basis and are used to define the convergent
Gaussian measure. The perturbative rules for the computation of expectation
values of graph-basis observables of any degree are given in terms of the
representation theoretic parameters. Explicit results for a number of
observables of degree up to four are given along with a Sage programme that
computes general expectation values.
