ARMA (p, q) model

An auto-regressive moving average model (ARMA(p, q)) combines the autoregression (AR(p)) and moving average (MA(q)) models into one. The current value of the simulated vector depends both on previous values of the same vector as well as previous values of the noise vector.

Complete the function definition of arma().

This exercise is part of the course

Optimizing R Code with Rcpp

View Course

Exercise instructions

Define an integer variable, start, equal to the maximum of p and q, plus one. Recall that max() is in the std namespace.
Inside the outer for loop, define a double variable, value, as mu plus the ith noise value.
Inside the first inner for loop, increase value by the jth element of theta times the "i minus j minus 1"th element of eps.
Inside the second inner for loop, increase value by the jth element of phi times the "i minus j minus 1"th element of x.

Hands-on interactive exercise

Have a go at this exercise by completing this sample code.

#include 
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector arma(int n, double mu, NumericVector phi, NumericVector theta, double sd) {
  int p = phi.size();
  int q = theta.size();   
  NumericVector x(n);
  
  // Generate the noise vector
  NumericVector eps = rnorm(n, 0.0, sd);
  
  // Start at the max of p and q plus 1
  ___
  
  // Loop i from start to n
  for(int i = start; i < n; i++) {
    // Value is mean plus noise
    ___
    
    // The MA(q) part
    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
      ___
    }
    
    // The AR(p) part
    for(int j = 0; j < p; j++) {
      // Increase by the jth element of phi times
      // the "i minus j minus 1"th element of x
      ___
    }
    
    x[i] = value;
  }
  return x;
}

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

Edit and Run Code