Here we demonstrate that tensor network techniques - originally devised for the analysis of quantum many-body problems - are well suited for the detailed study of rare event statistics in kinetically constrained models (KCMs). As concrete examples we consider the Fredrickson-Andersen and East models, two paradigmatic KCMs relevant to the modelling of glasses. We show how variational matrix product states allow to numerically approximate - systematically and with high accuracy - the leading eigenstates of the tilted dynamical generators which encode the large deviation statistics of the dynamics. Via this approach we can study system sizes beyond what is possible with other methods, allowing us to characterise in detail the finite size scaling of the trajectory-space phase transition of these models, the behaviour of spectral gaps, and the spatial structure and entanglement properties of dynamical phases. We discuss the broader implications of our results.