Mô phỏng mô hình MA(q)

Các mô hình moving average (MA) cũng phụ thuộc vào vòng lặp trước. Khác với mô hình AR, sự phụ thuộc nằm ở phần nhiễu.

Đây là thuật toán bằng 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 là số lượng quan sát được mô phỏng, mu là kỳ vọng, theta là một vector số của các hệ số moving average, và sd là độ lệch chuẩn của nhiễu.

Trước đó trong chương, bạn đã dùng R::rnorm() để sinh một số duy nhất từ phân phối chuẩn. Cũng có Rcpp::rnorm(), có thể sinh cả một vector số trong một lần gọi. Hàm này nhận cùng các đối số như rnorm() của R. Hoàn thiện phần định nghĩa hàm ma2(), bản dịch C++ của ma1().

Bài tập này là một phần của khóa học

Tối ưu hóa mã R với Rcpp

Xem khóa học

Hướng dẫn bài tập

Sinh vector nhiễu và gán vào eps. Dùng rnorm() từ namespace Rcpp (không phải namespace R).
Bên trong vòng lặp for phía ngoài, tính value bằng mu cộng với giá trị nhiễu thứ i.
Bên trong vòng lặp for phía trong, tăng value thêm phần tử thứ j của theta nhân với phần tử thứ "i trừ j trừ 1" của eps.
Sau các vòng lặp, gán phần tử thứ i của x bằng value.

Bài tập tương tác thực hành trực tiếp

Hãy thử làm bài tập này bằng cách hoàn thành đoạn mã mẫu này.

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

Chỉnh sửa và Chạy Mã