Simulate MA(q) model

Moving average (MA) models also depend upon the previous iteration. Unlike AR models, the dependency is on the noise part.

Here's the algorithm in 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 is the number of simulated observations, mu is the expected value, theta is a numeric vector of moving average coefficients, and sd is the standard deviation of the noise.

Earlier in the chapter, you used R::rnorm() to generate a single number from a normal distribution. There is also Rcpp::rnorm(), which can generate a whole numeric vector worth in one go. This takes the same arguments as R's rnorm(). Complete the function definition of ma2(), a C++ translation of ma1().

This exercise is part of the course

Optimizing R Code with Rcpp

View Course

Exercise instructions

Generate the noise vector as eps. Use rnorm() from the Rcpp namespace (not the R namespace).
Inside the outer for loop, calculate value as mu plus the ith noise value.
Inside the inner for loop, increase value by the jth element of theta times the "i minus j minus 1"th element of eps.
After the loops, set ith element of x to value.

Hands-on interactive exercise

Have a go at this exercise by completing this sample 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()
*/

Edit and Run Code