We have developed a new iterative geostatistical seismic amplitude variation with angle (AVA) inversion algorithm that inverts prestack seismic data, sorted by angle gathers, directly for high-resolution density, P-wave velocity, S-wave velocity, and facies models. This novel iterative geostatistical inverse procedure is based on two key main principles: the use of stochastic sequential simulation and cosimulation as the perturbation technique of the model parameter space and a global optimizer based on a crossover genetic algorithm to converge the simulated earth models toward an objective function, in this case, the mismatch between the recorded and synthetic prestack seismic data. As a geostatistical approach, all the elastic models simulated during the iterative procedure honors the well-log data at its own locations, the marginal prior distributions of P-wave velocity and S-wave velocity, and density estimated from the available well-log data, and the corresponding joint distributions between density versus P-wave velocity and P-wave versus S-wave velocity. We successfully tested and implemented this new algorithm on a synthetic prestack data set that mimicked the main properties of a real reservoir, and on a real seismic data set acquired over a deepwater turbidite oil reservoir. In both cases, the results showed a good convergence between the recorded and synthetic seismic. The synthetic example showed high-resolution inverted petroelastic models that reproduced the true petroelastic models. The inverted petroelastic models retrieved from the real case study found high resolution and do agree with previous seismic reservoir characterization studies.