  A multi--state life insurance model is naturally described in terms of the
intensity matrix of an underlying (time--inhomogeneous) Markov process which
describes the dynamics for the states of an insured person. Between and at
transitions, benefits and premiums are paid, defining a payment process, and
the technical reserve is defined as the present value of all future payments of
the contract. Classical methods for finding the reserve and higher order
moments involve the solution of certain differential equations (Thiele and
Hattendorf, respectively). In this paper we present an alternative
matrix--oriented approach based on general reward considerations for Markov
jump processes. The matrix approach provides a general framework for
effortlessly setting up general and even complex multi--state models, where
moments of all orders are then expressed explicitly in terms of so--called
product integrals (matrix--exponentials) of certain matrices. As Thiele and
Hattendorf type of theorems can be retrieved immediately from the matrix
formulae, this methods also provides a quick and transparent approach to
proving these classical results. Methods for obtaining distributions and
related properties of interest (e.g. quantiles or survival functions) of the
future payments are presented from both a theoretical and practical point of
view (via Laplace transforms and methods involving orthogonal polynomials).
