In this paper, we present some analytic properties of hidden variable bivariable fractal interpolation functions (HVBFIFs) with four function contractivity factors presented in [C. H. Yun and M. K. Li, Hidden variable bivariate fractal interpolation functions and errors on perturbations of function vertical scaling factors, Asian-Eur. J. Math. (2017), doi:10.1142/s1793557119500219]. Since four contractivity factors of these HVBFIFs are all functions, the construction of these HVBFIFs has more flexibility and diversity in fitting and approximation of complicated surfaces in nature and irregular experimental data with less self-similarity than one whose four contractivity factors are all constants or only one factor is function. The smoothness and stability of HVBFIFs are needed to ensure the applicability of the HVBFIFs in many practical problems such as the simulation of the objects of the nature, data fitting, etc. We first obtain the results related to their smoothness in nine different cases and then prove that the HVBFIFs are stable to the small perturbations of the interpolation points.