In this article, we discuss a model order reduction method for multiple-input and multiple-output discrete-time bilinear control systems. Similar to the continuous-time case, we will show that a system can be characterised by a series of generalised transfer functions. This will be achieved by a multivariate Z-transform of kernels corresponding to an explicit solution formula for discrete-time systems. We will further address the problem of generalised tangential interpolation which naturally comes along with this approach. We will introduce a reasonable generalisation of the linear H-2-norm. Based on this concept, we discuss the choice of interpolation points. Furthermore, an efficient discretisation of continuous-time systems is provided. The performance of the proposed method is evaluated in some numerical examples and compared with the method of balanced truncation for bilinear systems.