AR(p) モデルをシミュレーションする

自己回帰（Auto-regressive; AR）モデルは、時系列に対する線形回帰の一種で、予測値が過去の時点の値に依存します。過去の値に依存するため、モデルは時間を一つずつ進めて計算する必要があります。つまり for ループの出番で、C++ が役立ちます！

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 はシミュレーションする観測数、c は定数、phi は自己相関係数の数値ベクトル、eps はノイズの標準偏差です。ar1() を C++ に移植した ar2() の定義を完成させてください。

この演習はコースの一部です

Rcpp で R コードを最適化する

コースを見る

演習の手順

Rcpp の R API を使って、平均 c、標準偏差 eps の正規乱数を生成してください。
内側の for ループは 0 から p まで反復させます。
内側のループ内で、value を「phi の j 番目の要素 × x の「i から j を引いてさらに 1 引いた位置」の要素」だけ増やしてください。

実践的なインタラクティブ演習

このサンプルコードを完成させて、この演習に挑戦してみましょう。

#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()
*/

コードを編集して実行