ARMA stochastic process

Consider a stationary multivariate stochastic process X: \Omega \times[0,T] \rightarrow \Rset^d of dimension d, where X_t : \Omega \rightarrow \Rset^d is the random variable at time t. Under some general conditions,X can be modeled by an ARMA(p,q) model defined at the time stamp t by:

(1)X_t + \mat{A}_{\, 1} \, X_{t-1} + \hdots + \mat{A}_{\, p} \, X_{t-p} = \varepsilon_{t}+ \mat{B}_ {\, 1} \, \varepsilon_{t-1}+ \hdots + \mat{B}_{\, q} \, \varepsilon_{t-q}

where the coefficients of the recurrence are matrix in \Rset^d and (\varepsilon_t)_tis white noise discretized on the same time grid as the process X.

The coefficients (\mat{A}_{\, 1} , \hdots, \mat{A}_{\, p} ) form the Auto Regressive (AR) part of the model, while the coefficients (\mat{B}_{\, 1} \, \hdots, \mat{B}_{\, q} ) the Moving Average (MA) part.

We introduce the homogeneous system associated to (1):

(2)X_t + \mat{A}_{\, 1} \, X_{t-1} + \hdots + \mat{A}_{\, p} \, X_{t-p} = 0

To get stationary solutions of (1), it is necessary to get its characteristic polynomial defined in (3):

(3)\Phi(\vect{r}) = \vect{r}^p + \sum_{i=1}^p a_i\vect{r}^{p-i}

Thus the solutions of (2) are of the form P(t)\vect{r}_i^t where the \vect{r}_i are the roots of the polynomials \Phi(\vect{r}) defined in (3) and P is a polynomials of degree the order of the root \vect{r}_i:

The processes P(t)\vect{r}_i^t decrease with time if and only if the modulus of all the components of the roots \vect{r}_i are less than 1:

(4)\forall i,j \in [1,d], \, |r_{ij}| <1

Once given the coefficients of the model ARMA(p,q), we evaluate the roots of the polynomials \Phi(\vect{r}) and checks the previous condition (4). The roots \vect{r}_i, are the eigenvalues of the matrix \mat{M} which writes in dimension d as:

(5)\mat{M} = \left( \begin{array}{cccccc} \mat{0}_{\, d} & \mat{1}_{\, d} & \mat{0}_{\, d} & \hdots & \mat{0}_{\, d} & \mat{0}_{\, d} \\ \mat{0}_{\, d} & \mat{0}_{\, d} & \mat{1}_{\, d} & \hdots & \mat{0}_{\, d} & \mat{0}_{\, d}\\ \hdots \\ \mat{0}_{\, d} & \mat{0}_{\, d} & \mat{0}_{\, d} & \hdots & \mat{0}_{\, d} & \mat{1}_{\, d}\\ - \mat{A}_{\, 1} & -\mat{A}_{\, 2} & - \mat{A}_{\, 3}& \hdots & -\mat{A}_{\, p-1}& -\mat{A}_{\, p} \end{array} \right)

and in dimension 1:

(6)\mat{M} = \left( \begin{array}{cccccc} 0 & 1 & 0 & \hdots & 0 & 0\\ 0 & 0 & 1 & \hdots & 0 & 0\\ \hdots \\ 0 & 0 & 0 & \hdots & 0 & 1\\ -\alpha_1 & -\alpha_2 & -\alpha_3 & \hdots & -\alpha_{p-1} & -\alpha_p \end{array} \right)

The matrix \mat{M} is known to be the companion matrix.

It is important to note that:

  • when asking for a realization of the stationary process modeled by ARMA(p,q), one has to obtain a realization that does not depend on the current state of the process;

  • whereas, when one asks for a possible future extending a particular current state of the process, the realization of the model must depend on that current sate.

How to proceed to respect these constraints?

If we note \vect{X}_1(\omega,t) and \vect{X}_2(\omega,t) two distinct solutions of (1) associated to two distinct initial states, then, the process \vect{D}(\omega,t) = \vect{X}_2(\omega,t) - \vect{X}_1(\omega,t) is solution of the homogeneous equation associated to (1) and then decreases with time under the condition (4). Let us note N_{ther} the number such that:

(7)\left( \max_{i,j} |r_{ij}| \right)^{N_{ther}} < \varepsilon ,\Longleftrightarrow N_{ther} > \displaystyle \frac{\ln \varepsilon}{\ln \max_{i,j} |r_{ij}|}

where the r_i are the roots of the polynomials (3) and \varepsilon is the precision of the computer ( \varepsilon =2^{-53} \simeq 10^{-16}). Then, after N_{ther} instants, the process \vect{D}(\omega,t) has disappeared, which means that the processes \vect{X}_1(\omega,t) and \vect{X}_2(\omega,t) do not differ any more. As a conclusion, after N_{ther} instants, the realization of the ARMA process does not depend any more on the initial state.

That is why, when making a realization of the ARMA model, we perform a thermalization step that simply consists in realizing the model upon N_{ther} additional instants, erasing the N_{ther} first values and finally only retaining the other ones. That step ensures that the realization of the process does not depend on the initial state.

By default, the number N_{ther} is evaluated according to (7) by the method computeNThermalization. The User could get access to it with the method getNThermalization and can change the value with the method setNThermalization. (In order to give back to N_{ther} its default value, it is necessary to re-use the method computeNThermalization).

On the contrary, in the context of getting a possible future from a specified current state, the User should care that the number of additional instants N_{it} on which he wants to extend the process, is such that N_{it} < N_{ther} because beyond N_{ther}, the future has no link with the present. More precisely, after N_{it}^* instants, such that:

(8)\left( \max_{i,j} |r_{ij}| \right)^{N_{it}^*} < \max_{i} \sigma_i, \Longleftrightarrow N_{ther} > \displaystyle \frac{\max_{i} \sigma_i}{\ln \max_{i,j} |r_{ij}|}

where the \sigma_i are the components of the covariance matrix of the white noise \vect{\varepsilon}, the influence of the initial state is of same order than the influence of the white noise.

Let us note that when the ARMA model is created without specifying the current state, we automatically proceed to a thermalization step at the creation of the ARMA object.

Before asking for the generation of a possible future, the user has to specify the current state of the ARMA model, thanks to the creation method that takes into account the current state. In that case, we do not proceed to the thermalization step.

As an ARMA model is a stochastic process, the object ARMA inherits the methods of the Process object. Thus, it is possible to get its marginal processes, its time grid, its dimension and to get several realizations at a time of the process.