I don't have a reference, but here's the next best thing: a proof.
First, let's fix some notation. Let $P$ be a p-group, $G$ it's group of automorphisms, and $\Phi(P) = P^p[P,P]$ it's Frattini subgroup. Define inductively $\Phi^k(P)$ as $\Phi(\Phi^{k-1}(P)).$
The subgroups $\Phi^k(P)$ form a decreasing chain of characteristic subgroups which exhaust $P.$
Let $$G_k := \ker(G \rightarrow Aut(P/\Phi^k(P)))$$ and $$G_k' := \ker(G \rightarrow Aut(\Phi^k(P)/\Phi^{k+1}(P))).$$
Then the subgroups $G_k$ form an decreasing chain of normal subgroups of $G$ which exhaust $G.$
Let $d_k = \dim_{\mathbb{F}_p}(\Phi^k(P)/\Phi^{k+1}(P)).$ The group $P$ can be generated by $d_0$ elements. Choose a generating set $g_1 ... g_{d_0}$ and consider the map from $G_k \cap G_k'/G_{k+1}$ to $(\Phi^k(P)/\Phi^{k+1}(P))^{d_0}$
given by
$$\sigma \mapsto (\sigma(g_i)g_i^{-1})_{i=1}^{d_0}.$$
This map is an injective group homomorphism.
On the other hand $G_k/(G_{k} \cap G_{k}')$ injects into $Aut(\Phi^k(P)/\Phi^{k+1}(P)) \cong GL_{d_i}( \mathbb{F}_p).$
Note that if $p^n = |P|$ and $r = max\{d_i : d_i \neq 0\},$ then $\displaystyle\sum_{i=0}^r d_i = n.$ It follows that the order $$|G| = \displaystyle\prod_{k=0}^{r}|G_k/G_{k+1}| = \displaystyle\prod_{k=0}^{r} |(G_k \cap G_k')/G_{k+1}||(G_k \cap G_k')/G_{k+1}|$$ divides
$$\displaystyle\prod_{s=0}^{d_0-1} (p^{d_0} - p^s)\displaystyle\prod_{k=1}^{r} p^{d_kd_1}\displaystyle\prod_{s=0}^{d_k-1} (p^{d_k} - p^s)$$$$\displaystyle\prod_{s=0}^{d_0-1} (p^{d_0} - p^s)\displaystyle\prod_{k=1}^{r} p^{d_kd_0}\displaystyle\prod_{s=0}^{d_k-1} (p^{d_k} - p^s)$$
which divides
$$\displaystyle\prod_{k=0}^{n-1} (p^n - p^k).$$
