Obtaining the distribution of the profit and loss (PL) of a portfolio is a key problem in market risk measurement. However, existing methods, such as those based on the Normal distribution, and historical simulation methods, which use empirical distribution of risk factors, face difficulties in dealing with at least one of the following three problems: describing the distributional properties of risk factors appropriately (description problem); deriving distributions of risk factors with time horizon longer than one day (time aggregation problem); and deriving the distribution of the PL given the distributional properties of the risk factors (risk aggregation problem).Here, we show that expansion methods can provide reasonable solutions to all three problems. Expansion methods approximate a probability density function by a sum of orthogonal polynomials multiplied by an associated weight function. One of the most important advantages of expansion methods is that they only require moments of the target distribution up to some order to obtain an approximation. Therefore they have the potential to be applied in a wide range of situations, including in attempts to solve the three problems listed above. On the other hand, it is also known that expansions lack robustness: they often exhibit unignorable negative density and their approximation quality can be extremely poor. This limits applications of expansion methods in existing studies.In this thesis, we firstly develop techniques to provide robustness, with which expansion methods result in a practical approximation quality in a wider range of examples than investigated to date. Specifically, we investigate three techniques: standardisation, use of Laguerre expansion and optimisation. Standardisation applies expansion methods to a variable which is transformed so that its first and second moments are the same as those of the weight function. Use of Laguerre expansions applies those expansions to a risk factor so that heavy tails can be captured better. Optimisation considers expansions with coefficients of polynomials optimised so that the difference between the approximation and the target distribution is minimised with respect to mean integrated squared error. We show, by numerical examples using data sets of stock index returns and log differences of implied volatility, and GARCH models, that expansions with our techniques are more robust than conventional expansion methods. As such, marginal distributions of risk factors can be approximated by expansion methods. This solves a part of the description problem: the information on the marginal distributions of risk factors can be summarised by their moments. Then we show that the dependence structure among risk factors can be summarised in terms of their cross-moments. This solves the other part of the description problem. We also use the fact that moments of risk factors can be aggregated using their moments and cross-moments, to show that expansion methods can be applied to both the time and risk aggregation problems. Furthermore, we introduce expansion methods for multivariate distributions, which can also be used to approximate conditional expectations and copula densities by rational functions.