We study estimation of the parametric components of single and multiple index volatility models. Using the first- and second-order Stein’s identity, we develop methods that are applicable for estimation of the variance index in a high-dimensional setting requiring finite moment condition, which allows for heavy-tailed data. Our approach complements the existing literature in a low-dimensional setting, while relaxing the conditions on estimation, and provides a novel approach in a high-dimensional setting. We prove that the statistical rate of convergence of our variance index estimators consists of a parametric rate and a nonparametric rate, where the latter appears from the estimation of the mean link function. However, under standard assumptions, the parametric rate dominates the rate of convergence and our results match the minimax optimal rate for the mean index estimation. Simulation results illustrate finite sample properties of our methodology and back our theoretical conclusions.