Reachability analysis is a powerful tool for computing the set of states or outputs reachable for a system. While previous work has focused on systems described by state-space models, we present the first methods to compute reachable sets of ARMAX models - one of the most common input-output models originating from data-driven system identification. The first approach we propose can only be used with dependency-preserving set representations such as symbolic zonotopes, while the second one is valid for arbitrary set representations but relies on a reformulation of the ARMAX model. By analyzing the computational complexities, we show that both approaches scale quadratically with respect to the time horizon of the reachability problem when using symbolic zonotopes. To reduce the computational complexity, we propose a third approach that scales linearly with respect to the time horizon when using set representations that are closed under Minkowski addition and linear transformation and that satisfy that the computational complexity of the Minkowski sum is independent of the representation size of the operands. Our numerical experiments demonstrate that the reachable sets of ARMAX models are tighter than the reachable sets of equivalent state space models in case of unknown initial states. Therefore, this methodology has the potential to significantly reduce the conservatism of various verification techniques.