ARMA模型的引數估計方法

             ARMA引數估計和前面我們介紹的點估計內容相似，也介紹矩估計與最小二乘估計兩種方法。

           和上一次的點估計一樣，這一次我分享的內容主要有：矩估計，最小二乘估計，一個應用例題

            關於矩估計與最小二乘估計的基本思想，參見前面點估計的有關介紹.

           A RMA 模型（Auto-Regressive and Moving Average Model）是研究時間序列的重要方法，由自迴歸模型（簡稱AR模型）與滑動平均模型（簡稱MA模型）為基礎“混合”構成。在市場研究中常用於長期追蹤資料的研究，如：Panel研究中，用於消費行為模式變遷研究；在零售研究中，用於具有季節變動特徵的銷售量、市場規模的預測等。

           ARMA模型的基本原理：將預測指標隨時間推移而形成的資料序列看作是一個隨機序列，這組隨機變數所具有的依存關係體現著原始資料在時間上的延續性。一方面，影響因素的影響，另一方面，又有自身變動規律，假定影響因素為x1，x2，…，xk，由迴歸分析， 　　 

          我在這裡只介紹兩個根據資料直接估計出ar,arma,模型的函式:ar()，arima().其中ar(),arima()是擴充套件包TSA中的函式，如果要呼叫，你需要安裝並載入這個擴充套件包。

先看ar()的用法：

ar(x, aic = TRUE, order.max = NULL,

  method=c("yule-walker", "burg", "ols","mle", "yw"),

  na.action, series, ...)

這裡如果是對ar模型做矩估計，method選擇“yw”，極大似然估計使用”mle“。

我們來看一個例子：

還是用我們在上一篇博文中提到的生成ar模型資料的辦法生成一批資料，然後對他們做矩估計。

>w<-rnorm(550)

>x<-filter(w,filter=c(1,-0.9),"recursive")

>ar(x,order=2,method="yw")

Call:

ar(x = x,order.max = 2, method = "yw")

Coefficients:

      1       2 

 0.9870 -0.8905 

Orderselected 2  sigma^2 estimated as  1.085

對比係數，發現這個估計是十分接近的。這是因為對ar模型的矩估計利用的yule-waler方程可以看做線性方程，從而估計效果十分良好。可以對比最小二乘估計：

> ar(x,order=2,method="ols")

Call:

ar(x = x, order.max = 2, method ="ols")

Coefficients:

     1        2 

 0.9953 -0.8988 

Intercept: -0.003315 (0.04332)

Order selected 2  sigma^2 estimated as  1.028

可以看出估計效果差不多，同樣的，可以選擇極大似然估計，這裡不再贅述。

對於ma模型，必須指出的是無論是最小二乘，極大似然估計，還是矩估計，效果都是不理想的。在R中除了可以令arma(0,q)表示ma模型外（arma，arima提供的方法都是比較精細的，像矩估計這樣粗糙的方法是沒被採納的），我也沒有找到ma模型的矩估計函式（可能是因為效果不好。所以要用只有自己編寫的了）。

對於arma模型，可以使用函式arima(),arma().去做引數估計。這裡將兩個函式的用法摘錄如下：

arma(x, order = c(1, 1), lag = NULL, coef = NULL,

     include.intercept = TRUE, series = NULL, qr.tol = 1e-07, ...)

arima(x, order = c(0, 0, 0),

      seasonal = list(order = c(0, 0, 0), period = NA),

      xreg = NULL, include.mean = TRUE,

      transform.pars = TRUE,

      fixed = NULL, init = NULL,

      method = c("CSS-ML", "ML", "CSS"),

      n.cond, optim.method = "BFGS",

      optim.control = list(), kappa = 1e6)

最後舉一例來說明arma的引數估計：

還是用太陽黑子的資料來說明問題：

>data(sunspots)

>arma(sunspots,order=c(2,1))

Call:

arma(x =sunspots, order = c(2, 1))

Coefficient(s):

      ar1       ar2        ma1  intercept 

   1.1984   -0.2116    -0.6214     0.6673