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 increasing 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_1$ elements. Choose a generating set $g_1 ... g_{d_1}$ and consider the map
$$(G_{k} \cap G'\_{k})/G_{k+1} \rightarrow (\Phi^k(P)/\Phi^{k+1}(P))^{d_1}$$
given by
$$\sigma \mapsto (\sigma(g_i)g_i^{-1})_{i=1}^{d_1}.$$
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=1}^r d_i = n.$ It follows that the order $$|G| = \prod_{k=0}^{r}|G_k/G_{k+1}| = \prod_{k=0}^{r} |(G_k \cap G'\_k)/G_{k+1}||(G_k \cap G'\_k)/G_{k+1}|$$ divides
$$\displaystyle\prod_{k=0}^{r} p^{d_kd_1}\displaystyle\prod_{s=1}^{d_k-1} (p^{d_k} - p^s)$$
which divides
$$\displaystyle\prod_{k=0}^{n-1} (p^n - p^k).$$
