Value at Risk (VaR) has been successfully estimated using single covariate probabilistic fuzzy systems (PFS), a method which combines a linguistic description of the system behaviour with statistical properties of data. In this paper, we consider VaR estimation based on a PFS model for density forecast of a continuous response variable conditional on a high-dimensional set of covariates. The PFS model parameters are estimated by a novel two-step process. The performance of the proposed model is compared to the performance of a GARCH model for VaR estimation of the S&P 500 index. Furthermore, the additional information and process understanding provided by the different interpretations of the PFS models are illustrated. Our findings show that the validity of GARCH models are sometimes rejected, while those of PFS models of VaR are never rejected. Additionally, the PFS model captures both instant and periods of high volatility, and leads to less conservative models.