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.