This paper deals with non-Lipschitz nonlinear systems. Such systems can be approximated by a linear map of so-called signatures, which play a crucial role in the theory of rough paths and can be interpreted as collections of iterated integrals involving the control process. As a consequence, we identify a universal bilinear system, solved by the signature, that can approximate the state or output of the original nonlinear dynamics arbitrarily well. In contrast to other (bi)linearization techniques, the signature approach remains feasible in large-scale settings, as the dimension of the associated bilinear system grows only with the number of inputs. However, the signature model is typically of high order, requiring an optimization process based on model order reduction (MOR). We derive an MOR method for unstable bilinear systems with non-zero initial states and apply it to the signature, yielding a potentially low-dimensional bilinear model. An advantage of our method is that the original nonlinear system need not be known explicitly, since only data are required to learn the linear map of the signature. The subsequent MOR procedure is model-oriented and specifically designed for the signature process. Consequently, this work has two main applications: (1) efficient modeling/data fitting using small-scale bilinear systems, and (2) MOR for nonlinear systems. We illustrate the effectiveness of our approach in the second application through numerical experiments.