This study uses Genetic Programming (GP) to discover new types of volatility forecasting models for financial time series. GP is a convenient tool to explore the space of potential forecasting models and to select the more robust solutions. The application to foreign exchange financial problems requires an exact symmetry induced by the interchange of currencies. In GP, this symmetry is enforced by using a strongly typed GP approach and syntactic restrictions on the node set. GP convergence is increased by a few orders of magnitude by optimizing the constants in the GP trees with a local optimization algorithm. The various algorithms are compared on the discovery of the symmetric transcendental function cosine. For the volatility forecast, the optimization is performed using return time series sampled hourly, possibly including aggregated returns at longer time horizons. The in-sample optimization and out-of-sample tests are performed on 13 years of high frequency data for two foreign exchange time series. The out-of-sample forecasting performance of these new models are compared with the corresponding performance of some popular ARCH-types models, and GP consistently outperform the benchmarks. In particular, GP discovered that cross products of returns at different time horizons improve substantially the forecasting performance.