Recent rapid advances in single particle tracking and supercomputing techniques resulted in an unprecedented abundance of diffusion data exhibiting complex behaviours, such the presence of power law tails of the msd and memory functions, commonly referred to as "fractional dynamics". Motivated by these developments, we study the stationary solutions of the classical and generalized Langevin equation as models of the contemporarily observed phenomena. In the first chapter we sketch the historical background of the generalized Langevin equation. The second chapter is devoted to the brief overview of the theory of the Gaussian variables and processes. In the third chapter we derive the generalized Langevin equation from Hamilton's equations of motion. In the fourth chapter the series of propositions and theorems shows that a large class of Langevin equations has a solutions that, sampled in discrete time, are the moving average autoregressive processes, which can be analysed using large number of available statistical methods. The fifth chapter starts with a short introduction into the basic notions of ergodic theory. Then we prove a generalized version of the classical Maruyama's theorem, which describes the time averages of all stationary Gaussian processes with a finite number of spectral atoms. In the last, sixth chapter we continue the study of the non-ergodic solutions of the Langevin equation, this time introducing the ergodicity breaking through the random parametrisation of the equation itself. This leads to a specific form of non-Gaussianity which affects only many-dimensional distributions of the process. The main result of this chapter is the series of propositions which describes the second order structure of solutions of the superstatistical Langevin equation.