Computing statistical quantities of interest of the solution of PDE on random domains is an important and challenging task in engineering. We consider the computation of these quantities by the perturbation approach. Especially, we discuss how third order accurate expansions of the mean and the correlation can numerically be computed. These expansions become even fourth order accurate for certain types of boundary variations. The correction terms are given by the solution of correlation equations in the tensor product domain, which can efficiently be computed by means of H-matrices. They have recently been shown to be an efficient tool to solve correlation equations with rough data correlations, that is, with low Sobolev smoothness or small correlation length, in almost linear time. Numerical experiments in three dimensions for higher order ansatz spaces show the feasibility of the proposed algorithm. The application to a non-smooth domain is also included.