  We present the results of an exploratory study of the numerical stochastic
perturbation theory (NSPT) applied to the four dimensional twisted Eguchi-Kawai
(TEK) model. We employ a Kramers type algorithm based on the Generalized Hybrid
Molecular Dynamics (GHMD) algorithm. We have computed the perturbative
expansion of square Wilson loops up to $O(g^8)$. The results of the first two
coefficients (up to $O(g^4)$) have a high precision and match well with the
exact values. The next two coefficients can be determined and even extrapolated
to large $N$, where they should coincide with the corresponding coefficients
for ordinary Yang-Mills theory on an infinite lattice. Our analysis shows the
behaviour of the probability distribution for each coefficient tending to
Gaussian for larger $N$. The results allow us to establish the requirements to
extend this analysis to much higher order.
