Wheat is one of the most important global crops, understanding the drivers of wheat yield has significant societal benefits. Climate variables are particularly important in determining interannual variations in wheat yield, either as primary factors which directly influence the stages of wheat growth, or as secondary factors through their influence on pests, diseases and soil conditions. Here we present a new approach to model wheat yield; an empirical method based on nonlinear complex systems identification, known as NARMAX (Nonlinear AutoRegressive Moving Average with eXogenous inputs model). We deploy the NARMAX analytical approach for a specific site, Rothamsted, UK, where detailed meteorological variables are available, together with specific information on site conditions and crop growth stages. NARMAX yield forecasts are compared with those from the WOFOST crop model and nine state-of-the-art machine learning (ML) models; experimental results show that NARMAX outperforms all the compared methods in both prediction accuracy and model interpretability. We also develop regional wheat yield forecasts derived from a new gridded meteorological data product.The NARMAX approach produces skillful forecasts (r = 0.78) of Rothamsted wheat yield for a validation period, with small errors. The NARMAX regional forecasts, based on less specific information than WOFOST, also show a high degree of skill (r = 0.73). In addition, the predictor terms chosen for the model are identifiable and can help to give insight into potential key processes involved in the determination of wheat yield at a specific location. This approach can be extended in principle to other crop types and locations. It is straightforward and inexpensive to implement, using a limited number of meteorological predictor variables, which can be taken from site-based observations, or from gridded meteorological datasets. The method is a new tool to understand the environmental drivers of wheat yields on an annual basis.