Probit models with spatial dependencies were first studied by McMillen (1992), where an EM Algorithm was developed to produce consistent (maximum likelihood) estimates for these models. In spatial autoregressive probit model, the spatial dependent structure adds complexity in the estimation of parameters. LeSage and Smith (2001) use Bayesian estimation via Markov Chain Monte Carlo methods that sample sequentially from the complete set of conditional distributions for all parameters. Klier and McMillen (2008) have proposed a linearized version of the GMM estimator that avoids the infeasible problem of inverting n-by-n matrices when employing large samples. They show that standard GMM reduces to a nonlinear two-stage least squares problem. Martinetti and Geniaux (2017) proposed approximate likelihood estimation which based on the full maximization of likelihood of an approximate multivariate normal distribution function. We use some extensive simulation for these methods and show the best estimation method which can handle sample sizes with many observations and various value of coefficient spatial lag, provided the spatial weight matrix is inconvenient sparse form, as is for large data sets, where each observation neighbors only a few other observations.