In this paper, we discuss probit model on multivariate binary response. We assume that each of n individuals is observed in T responses. Yit is tth response on ith individual/subject and each response is binary. Each subject has covariate Xi (individual characteristic) and covariate Zijt (characteristic of alternative j). Response on individual ith can be represented by Yi = (Yi1,....,YiT), Yit is tth response on ith individual/subject and each response is multinomial. In order to simplify, we choose one of individual characteristics and alternative characteristics. We use simulated maximum likelihood estimator (SMLE) methods to estimate the parameter based on Geweke-Hajivassiliou-Keane (GHK) simulator. We find the first derivative of likelihood function for multivariate binary probit. Then, we expand to multivariate multinomial response. The first derivative is used in the BHHH (Berndt, Hall, Hall, Hausman) iteration to obtain estimators.