Simulate AR(p) model

Auto-regressive (AR) models are a type of linear regression for time series where the predicted values depend upon the values at previous time points. Due to this dependence on previous time points, the model must be calculated one time point after another. That means using a for loop, so C++ comes in handy!

Here's the algorithm in R:

ar1 <- function(n, constant, phi, eps) {
  p <- length(phi)
  x <- numeric(n)
  for(i in seq(p + 1, n)) {
    value <- rnorm(1, constant, eps)
    for(j in seq_len(p)) {
      value <- value + phi[j] * x[i - j]
    }
    x[i] <- value
  }
  x
}

n is the number of simulated observations, c is a constant, phi is a numeric vector of autocorrelation coefficients, and eps is the standard deviation of the noise. Complete the definition of ar2(), a C++ translation of ar1().

This exercise is part of the course

Optimizing R Code with Rcpp

View Course

Exercise instructions

Generate a normally distributed random number with mean c and standard deviation eps, using Rcpp's R API.
Make the inner for loop iterate from 0 to p.
Inside the inner 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 ar2(int n, double c, NumericVector phi, double eps) {
  int p = phi.size();  
  NumericVector x(n);
  
  // Loop from p to n
  for(int i = p; i < n; i++) {
    // Generate a random number from the normal distribution
    double value = ___::___(___, ___);
    // Loop from zero to p
    for(int j = ___; j < ___; j++) {
      // Increase by the jth element of phi times 
      // the "i minus j minus 1"th element of x
      value += ___[___] * ___[___];
    }
    x[i] = value;
  }
  return x;
}

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

Edit and Run Code