This thesis considers matrix methods in multi-state life insurance, with an emphasis on techniques related to inhomogeneous phase-type distributions (IPH) andproduct integrals. We start out with developing an expectation-maximization (EM)algorithm for statistical estimation of general IPHs. Then we introduce a newclass of multi-state models, the so-called aggregate Markov model, which allowsfor non-Markovian modeling with most of the analytical tractability of Markovchains preserved. Using techniques related to IPHs, we derive distributional properties, computational schemes for life insurance valuations with duration-dependentpayments, and statistical estimation procedures based on the EM algorithm forgeneral IPHs. Special attention is given to a case with a reset property, wherethe aggregate Markov model is semi-Markovian. We then move on and considerMarkov chain interest rate models and show that bond prices are survival functionsof IPHs. This allows for calibration via EM algorithms for phase-type distributions.Then we consider a multivariate payment process and derive higher order momentsof its present value. Finally, we consider computation of market values of bonuspayments in multi-state with-profit life insurance, where numerical procedures basedon simulation of financial scenarios and classic analytical methods for insurancerisk are developed.