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 est la valeur attendue, theta est un vecteur numérique de coefficients de moyenne mobile et sd est 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 distribution normale. Il existe aussi Rcpp::rnorm(), qui peut générer d'un coup un vecteur numérique complet. Cette fonction prend les mêmes arguments que rnorm() de R. Complétez la définition de la fonction ma2(), une traduction C++ de ma1().

Cette activité fait partie du cours

Optimiser le code R avec Rcpp

Voir le cours

Instructions de l’exercice

Générez le vecteur de bruit sous le nom 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 valeur de bruit en position i.
Dans la boucle interne, augmentez value du je élément de theta multiplié par l'élément « i moins j moins 1 » de eps.
Après les boucles, assignez l'élément en position i de x à value.

Exercice interactif pratique

Essayez cet exercice en complétant ce code d’exemple.

#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