Neural networks have been widely used to solve complex real-world problems. Due to the complicate, nonlinear, non-convex nature of neural networks, formal safety guarantees for the behaviors of neural network systems will be crucial for their applications in safety-critical systems. In this paper, the reachable set estimation and verification problems for Nonlinear Autoregressive-Moving Average (NARMA) models in the forms of neural networks are addressed. The neural network involved in the model is a class of feed-forward neural networks called Multi-Layer Perceptron (MLP). By partitioning the input set of an MLP into a finite number of cells, a layer-by-layer computation algorithm is developed for reachable set estimation for each individual cell. The union of estimated reachable sets of all cells forms an over-approximation of reachable set of the MLP. Furthermore, an iterative reachable set estimation algorithm based on reachable set estimation for MLPs is developed for NARMA models. The safety verification can be performed by checking the existence of intersections of unsafe regions and estimated reachable set. Several numerical examples are provided to illustrate our approach.