This paper presents a homotopy approach to improving PEM identification of ARMAX model. PEM estimates of ARMAX model parameters are determined as the global minimum of criterion function, which is however not always unimodal because of the MA noise model part. An optimization-based PEM identification algorithm has a potential risk to be stuck at a local minimum that results in a poorly identified model. A homotopy continuation method is introduced to solve this problem. The idea is to start the estimation with the criterion function for PEM identification of the ARX model, which is gradually deformed into the actual one for PEM identification of the ARMAX model as the algorithm iterates. By building the deformation into the usual recursive procedure for the ARMAX identification and introducing a scheme to control the solution continuously staying in the global minima of the deformed criterion functions, the homotopy-based PEM identification algorithm is implemented in such a way that it has very good convergence performance, with only little increase in computation load compared to the usual PEM algorithm.