Simuler un modèle MA(q)

Les modèles à moyenne mobile (MA) dépendent eux aussi de l’itération précédente. Contrairement aux modèles AR, la dépendance porte sur la partie bruit.

Voici l’algorithme en R :

ma1 <- function(n, mu, theta, sd) {
  q <- length(theta)
  x <- numeric(n)
  eps <- rnorm(n, 0, sd)
  for(i in seq(q + 1, n)) {
    value <- mu + eps[i]
    for(j in seq_len(q)) {
      value <- value + theta[j] * eps[i - j]
    }
    x[i] <- value
  }
  x
}

n est le nombre d’observations simulées, mu la valeur attendue, theta un vecteur numérique de coefficients de moyenne mobile, et sd l’écart type du bruit.

Plus tôt dans le chapitre, vous avez utilisé R::rnorm() pour générer un seul nombre selon une loi normale. Il existe aussi Rcpp::rnorm(), qui peut générer en une fois tout un vecteur numérique. Elle prend les mêmes arguments que rnorm() de R. Complétez la définition de la fonction ma2(), une traduction C++ de ma1().

Cet exercice fait partie du cours

Optimiser du code R avec Rcpp

Afficher le cours

Instructions

Générez le vecteur de bruit eps. Utilisez rnorm() de l’espace de noms Rcpp (et non de l’espace de noms R).
Dans la boucle externe, calculez value comme mu plus la i-ième valeur de bruit.
Dans la boucle interne, augmentez value du j-ième élément de theta multiplié par l’élément « i moins j moins 1 » de eps.
Après les boucles, affectez value au i-ième élément de x.

Exercice interactif pratique

Essayez cet exercice en complétant cet exemple de code.

#include 
using namespace Rcpp ;

// [[Rcpp::export]]
NumericVector ma2( int n, double mu, NumericVector theta, double sd ){
  int q = theta.size(); 
  NumericVector x(n);
  
  // Generate the noise vector
  NumericVector eps = ___(___, 0.0, ___);
    
  // Loop from q to n
  for(int i = q; i < n; i++) {
    // Value is mean plus noise
    double value = ___ + ___;
    // Loop from zero to q
    for(int j = 0; j < q; j++) {
      // Increase by the jth element of theta times
      // the "i minus j minus 1"th element of eps
      value += ___ * ___;
    }
    // Set ith element of x to value
    ___ = ___;
  }
    return x ;
}

/*** R
d <- data.frame(
  x = 1:50,
  y = ma2(50, 10, c(1, -0.5), 1)
)
ggplot(d, aes(x, y)) + geom_line()
*/

Modifier et exécuter le code