I want a procedure for randomly generating a square matrix A such that the linear system
x_t+1 = A x_t
is globally asymptotically stable (e.g. the discrete Lyapunov equation is satisfied) for a state vector x_t of arbitrary dimensionality.
The statistical properties of the matrices are not particularly important (though I am imagining that using a zero-mean Gaussian distribution somewhere in this process may be helpful).
Can somebody help?
A linear system is strongly stable if all its eigenvalues lie inside the unit disk. Presumably you want all entries in the linear system to be real, so any complex eigenvalues should come in complex-conjugate pairs.
If no other constraints are imposed you can simply reverse engineer the whole process by first selecting those eigenvalues randomly inside the unit disk, then building your matrix by either strategy below:
(1) select a random collection of column vectors to make a square matrix $P$ which will be the eigenvectors, and setting $ A= PDP^{-1}$. Here $D$ is the diagonal matrix that contains the eigenvalues.
If you want to avoid having to compute $P^{-1}$ and have no need to know what are the eigenvectors in $P$ you could instead
(2) create the companion matrix $A$ associated to the characteristic polynomial whose roots were chosen above. This matrix will have the desired eigenvalues. It will automatically be real, saving you from having to work with complex numbers.
Of course you can also perform a random change of basis on the matrix constructed in (2) if you want to obscure its structure. That means computing $BAB^{-1}$ where $B$ is any invertible matrix you choose to select.
Here is a toy example illustrating the steps in (2).