ARMA(1,1)和AR(2)自相关函数的R语言实现

ARMA(1,1)模型

ARMA(1,1)模型的自相关函数

  模型定义的方程是：

$$Y_t=\phi Y_{t-1}+e_t-\theta e_{t-1}$$

为得到Yule-Walk形式的方程组，首先：

$$E(e_tY_t)=E(e_t(\phi Y_{t-1}+e_t-\theta e_{t-1}))=\sigma^2_e$$

或

$$E(e_{t-1}Y_t)=E(e_{t-1}(\phi Y_{t-1}+e_t-\theta e_{t-1}))=(\phi-\theta)\sigma^2_e$$

在其定义式两边乘以$Y_{t-k}$并求期望值，得到：

$$\left.\begin{aligned} \gamma_0 & =\phi\gamma_1+[1-\theta(\phi-\theta)\sigma^2_e]\\ \gamma_1 & = \phi\gamma_0-\theta\sigma^2_e\\ \gamma_k & =\phi\gamma_{k-1} \end{aligned} \right \}$$

求解前两个方程得到：

$$\gamma_0=\frac{1-2\phi\theta+\theta^2}{1-\phi^2}\sigma^2_e$$

进一步通过递推关系式，得到:

$$\rho_k=\frac{(1-\theta\phi)(\phi-\theta)}{1-2\phi\theta+\theta^2}\phi^{k-1}，k>1$$

  值得注意的是，随着之后长度$k$的增加，该模型的相关函数指数递减。相关函数的初始值 $\rho_0=0$，从 $\rho_k$ 的递推式中可以看到，它的形状依赖于 $\rho_1$ 和 $\phi$ 的符号。

ARMA(1,1)模型自相关函数的R语言函数

  按照前述中的递推式可构造求解ARMA(1,1)模型自相关函数并画出其自相关图像的函数 ARMA1.1 <- function(phi=phi, theta=theta, max.lag=20, pic=T, ylim=NULL)

1
2
3
4
5
6
7
8
9
10
11
12
13

ARMA1.1 <- function(phi=phi, theta=theta, max.lag=20, pic=T, ylim=NULL){
  rho = NULL
  rho[1] <- 1
  for (k in 1:max.lag) {
    rho[k+1] <- (1-theta*phi)*(phi-theta)/(1-2*theta*phi+theta^2)*phi^(k-1)
  }
  if(pic == T){
    plot(x = 1:max.lag, y = rho[-1], type = "h", 
         xlab = "Lag", ylab = "ACF", ylim = ylim)
    abline(h=0)
  }
  out <- list(phi = phi, theta = theta, rho = rho)
}

  该函数包含五个参数:

1. phi  ARMA(1,1)模型中自回归部分的系数，即$\phi$.

2. theta  ARMA(1,1)模型中滑动平均部分的系数，即$\theta$.

3. max.lag  ARMA(1,1)模型中需要求解的最大滞后长度k，默认为$k=20$.

4. pic  逻辑值，pic=T时绘制自相关函数图像，默认pic=T.

5. ylim  一个二维的向量，作为绘制自相关函数图像的纵坐标范围，默认为ylim=NULL

该函数的输出为一个列表，其中phi为输入的 phi，theta 为输入的 theta，rho 为计算所得的max.lag+1个自相关函数（第一个为 $\rho_0=1$）。

AR(2)模型

AR(2)模型的自相关函数

  模型定义的方程是：

$$Y_t=\phi_1Y_{t-1}+\phi_2Y_{t-2}+e_t$$

为推导其自相关函数，在其两边乘上$Y_{t-k}$，并求期望。假定平稳、零均值，并且$e_t$独立于$Y_{t-1}$,得到：

$$\gamma_k=\phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}，k=1,2,3,\cdots$$

两边再除以$\gamma_0$，

$$\rho_k=\phi_1\rho_{k-1}+\phi_2\rho{k-2}，k=1,2,3,\cdots$$

称上述方程为Yule-Walker方程。

  下面需求解 $k=1$ 和 $k=2$ 时对应的$\rho$值。首先令 $k=1$，因 $\rho_0=1，\rho_{-1}=\rho_1$，得到 $\rho_1=\phi_1+\phi_2\rho_2$ ，整理得：

$$\rho_1=\frac{\phi_1}{1-\phi_2}$$

再令 $k=2$，将上式带入Yule-Walker方程中，得出

$$\rho_2=\frac{\phi_2(1-\phi_2)+\phi_1^2}{1-\phi_2}$$

至此我们可以通过Yule-Walker方程求出AR(2)模型的自相关函数。

  下面在给出计算$\rho_k$的显示表达：

$$\rho_k=\frac{(1-G^2_2)G^{k+1}_1-(1-G^2_1)G^{k+1}_2}{(G_1-G_2)(1+G_1G_2)}，k\ge 0$$

其中，$G1=\frac{\phi-\sqrt{\phi^2_1+4\phi_2}}{2}$、$G2=\frac{\phi+\sqrt{\phi^2_1+4\phi_2}}{2}$

特别地，当$\phi^2_1+4\phi_2<0$时，其可表示为：

$$\rho_k=R^k\frac{sin(\theta k+\Phi)}{sin(\phi)}，k\ge 0$$

其中，$R=\sqrt{-\phi_2}$ 为阻尼因子，$\theta=arccos(\phi_1/(2\sqrt{-\phi_2}))$ 为频率，$\phi=arctan[(1-\phi_2)/(1+\phi_2)]$

AR(2)模型自相关函数的R语言函数

  按照前述中的递推式可构造求解AR(2)模型自相关函数并画出其自相关图像的函数 AR2 <- function(phi1=phi1, phi2=phi2, max.lag = 20, pic = T, ylim=c(-1,1))，此函数还可求出其特征根为复数时的阻尼因子$R$和$\Theta$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

AR2 <- function(phi1=phi1, phi2=phi2, max.lag = 20, pic = T, ylim=c(-1,1)){
  rho = NULL;root=NULL;
  rho1 <- phi1/(1-phi2)
  rho2 <- (phi2*(1-phi2)+phi1^2)/(1-phi2)
  rho[1] <- rho1
  rho[2] <- rho2
  
  for (k in 3:max.lag) {
    rho[k] <- phi1*rho[k-1]+phi2*rho[k-2]
  }
  
  if(pic == T){
    plot(x = 1:max.lag, y = rho, type = "h", 
         xlab = "Lag", ylab = "ACF", ylim = ylim)
    abline(h=0)
  }
  
  if(phi1^2+4*phi2 < 0){
    root1 <- (phi1+sqrt(as.complex(phi1^2+4*phi2)))/(-2*phi2)
    root2 <- (phi1-sqrt(as.complex(phi1^2+4*phi2)))/(-2*phi2)
    R <- sqrt(-phi2)
    Theta <- acos(phi1/2*sqrt(-phi2))
  }else{
    root1 <- (phi1+sqrt(phi1^2+4*phi2))/(-2*phi2)
    root2 <- (phi1-sqrt(phi1^2+4*phi2))/(-2*phi2)
    R <- NULL
    Theta <-NULL
  }
  out <- list(phi=c(phi1,phi2), rho = rho, root = c(root1, root2), 
              R = R, Theta = Theta)
}

  该函数包含五个参数:

1. phi1,phi2  AR(2)模型中的两个参数，即$\phi_1，\phi_2$.

2. max.lag  AR(2)模型中需要求解的最大滞后长度k，默认为$k=20$.

3. pic  逻辑值，pic=T时绘制自相关函数图像，默认pic=T.

4. ylim  一个二维的向量，作为绘制自相关函数图像的纵坐标范围，默认为ylim=c(-1,1)

该函数的输出为一个列表，其中phi为输入的phi为输入的两个参数，rho为计算所得的max.lag个自相关函数，root为特征方程的两个根，R为阻尼系数，Theta为频率。