This paper proposes a three-stage procedure for parametric identification of piecewise affine autoregressive exogenous (PWARX) models. The first stage simultaneously classifies the data points and estimates the number of submodels and the corresponding parameters by solving the partition into a minimum number of feasible subsystems (MIN PIPS) problem for a suitable set of linear complementary inequalities derived from data. Second, a refinement procedure reduces misclassifications and improves parameter estimates. The third stage determines a polyhedral partition of the regressor set via two-class or multiclass linear separation techniques. As a main feature, the algorithm imposes that the identification error is bounded by a quantity delta. Such a bound is a useful tuning parameter to trade off between quality of fit and model complexity. The performance of the proposed PWA system identification procedure is demonstrated via numerical examples and on experimental data from an electronic component placement process in a pick-and-place machine.