In this paper, we propose an iterative generalized eigenvector (GE) precoder design for three dimensional (3D) massive multi-input multi-output (MIMO) downlink with uniform planar array (UPA) and imperfect channel state information (CSI). We use a posterior beam based statistical channel model which includes the channel aging and the spatial correlation. We consider the problem of maximizing the expected sum-rate under a total power constraint. By replacing the expected sum-rate with their deterministic equivalents, we obtain the structure of the optimal linear precoders. The columns of the optimal precoding matrices are the generalized eigen-vectors of two matrices. One of the two matrices is related to the transmitted signal of the intended user, and the other is related to the power constraint and the leakage to other users. Furthermore, the two matrices are also functions of the precoders, thus iterative updates of the precoder are needed. Then, we investigate the power allocation and the computation of the Lagrangian multiplier, and propose an iterative algorithm based on the structure of the optimal precoders. Simulation results show that the proposed precoders can achieve significantly performance gain than the widely used regularized zero forcing (RZF) precoder and signal to leakage plus noise ratio (SLNR) precoder.