Two methods of analytical approximations for computing the value of a European option on the conditional variance in a GARCH setting are presented. The first is based on the Johnson density functions family while the second uses the generalized Edgeworth expansion of the conditional variance’s risk-neutral density function. The analytical approximations based on the lognormal system of Johnson densities is found to be the most accurate one developed in this paper. However, the precision of the approximation is weaker for deep out-of-the-money options when the conditional variance dynamics of the NGARCH process displays a high level of persistence or almost no variability. When the moments of the conditional-variance dynamics display explosive behaviour, the precision of the approximation is generally good, except for the case of long-maturity options. The numerical analysis also suggests that the error term associated with the generalized Edgeworth expansion will have a critical impact when the skewness and kurtosis to be reproduced deviate, even slightly, from those of the approximate distribution used in the expansion. The proposed methodology can easily be adapted to other GARCH processes and generalized to the case of a volatility option.