自我回歸模型

用 Rcpp 最佳化 R 程式碼

Romain François

Consulting Datactive, ThinkR

自我回歸（AR）

ar <- function(n, phi, sd){
    x  <- epsilon <- rnorm(n, sd = sd)
    np <- length(phi)

    for( i in seq(np+1, n)){
        x[i] <- sum(x[seq(i-1, i-np)] * phi) + epsilon[i]
    }
    x
}

C++ 的 AR

先來第 1 步，填滿前 np 個值

NumericVector x(n) ;
// initial loop
for( ___ ; __ < np ; ___ ){
    x[i] = R::rnorm(___) ; }

主體包含外層與內層第 {{1}} 步

// outer loop
for( ___ ; ___ ; ___ ){    
    double value = rnorm(___) ;    
    // inner loop
    for( ___ ; ___ ; ___ ){
        value += ___ ;    }
    x[i] = value ;
}

移動平均模擬

ma <- function(n, theta, sd){
    epsilon <- rnorm(n, sd = sd)
    x <- numeric(n)
    nq <- length(theta)

    for( i in seq(nq+1, n)){
        x[i] <- sum(epsilon[seq(i-1, i-nq)] * theta) + epsilon[i]
    }
    x
}

移動平均

#include <Rcpp.h>
using namespace Rcpp ;
// [[Rcpp::export]]
NumericVector ma( int n, double mu, NumericVector theta, double sd ){
    int nq = theta.size() ;

    // 一次產生噪音向量
    // 使用 Rcpp::rnorm，類似 R 的函式
    NumericVector eps = Rcpp::rnorm(n, 0.0, sd) ;

    // 初始化長度為 n、值為 0.0 的輸出向量 
    NumericVector x(___) ;

    // 從索引 nq + 1 開始填值
    for( int i=nq+1; i<n; i++){
        ____   }
    return x ;
}

ARMA(p,q) = AR(p) + MA(q)

一起來練習吧！

用 Rcpp 最佳化 R 程式碼