Liouville Conformal Field Theory (LCFT) is an essential building block of Polyakov's formulation of non critical string theory. Moreover, scaling limits of statistical mechanics models on planar maps are believed by physicists to be described by LCFT. A rigorous probabilistic formulation of LCFT based on a path integral formulation was recently given by the present authors and F. David in \cite{DKRV}. In the present work, we prove the validity of the conformal Ward identities and the Belavin-Polyakov-Zamolodchikov (BPZ) differential equations (of order $2$) for the correlation functions of LCFT. This initiates the program started in the seminal work of Belavin-Polyakov-Zamolodchikov \cite{BPZ} in a probabilistic setup for a non-trivial Conformal Field Theory. We also prove several celebrated results on LCFT, in particular an explicit formula for the 4 point correlation functions (with insertion of a second order degenerate field) leading to a rigorous proof of a non trivial functional relation on the 3 point structure constants derived earlier in the physics literature by Teschner \cite{Tesc}. The proofs are based on exact identities which rely on the underlying Gaussian structure of LCFT combined with estimates from the theory of critical Gaussian Multiplicative Chaos and a careful analysis of singular integrals (Beurling transforms and generalizations). As a by-product, we give bounds on the correlation functions when two points collide making rigorous certain predictions from physics on the so-called "operator product expansion" of LCFT.