Parametric volatility models can be seen as the result of some form of dimensionality reduction obtained by projecting the volatility surface into a basis of risk factors. Examples include polynomial models and stochastic volatility models having an explicit expression for the smile, such as the SVI and the SABR model. When modelling the dynamics of the risk factors, one can define some stochastic processes, consider a statistical model, or assume that they are function of some exogenous explanatory variables (such as the spot price and volume). In the first two, we are tied to our ability of capturing the true dynamics of the underlying stochastic processes, leading to computation errors, in the latter we assume known the explanatory variables. Alternatively, we can use machine learning to implicitly infer the dynamics of the risk factors, or that of their explanatory variables, from market data. We let the risk factors be function of some stochastic explanatory variables: the spot price, the set of time-to-maturity, the trading volume, some external noise factors, the estimated volatility, and financial indicators like the VIX. We use a statistical model to capture the relation between the risk factors and their explanatory variables, and estimate empirically the model parameters. The risk factors being stochastic, we estimate the future implied volatility surface on average, by solving its conditional expectation with respect to the explanatory variables. This is a multi-step prediction problem, and we propose to use temporal difference backpropagation (TDBP) models for learning to predict the value function. This way, we use market data to estimate the implicit stochastic processes driving the dynamics of the future volatility surface. We test the performance of our model on three well known parametric volatility models against Monte Carlo simulations, and we show the effects of the explanatory variables on the dynamics of the volatility surface. In all cases, we recover nearly exactly the expected volatility surface generated by the Monte Carlo simulations with only one five hundredth paths. Our framework can be applied to option pricing and hedging as well as to risk management.