YY分享

1. ARMA模型

ARMA模型[1]是由美国统计学家G.E.P.Box和英国统计学家G.M.Jenkins在二十世纪七十年代提出的时序分析模型，即自回归滑动平均模型(Autoregressive Moving Average Model)，用此模型所作的时间序列预测方法也称博克斯-詹金斯(B-J)法。

2. 模型表达式

ARMA模型可以细分为自回归模型（auto regression model，AR模型）、滑动平均模型（moving average，MA模型）和ARMA模型三大类。

（1）p阶AR模型

（2）q阶MA模型

（3）ARMA模型

运用ARMA模型的前提条件是,建立模型的时间序列是由一个零均值的平稳随机过程产生的。即其过程的随机性质具有时间上的不变性,在图形上表现为所有的样本点皆在某一水平线上下随机地波动。

3. 实现算法介绍

其中ai (i=1, 2, …, p) 是自回归系数，bi (i=0, 1, …, q) 是滑动平均系数，w(n)是白噪声。给定白噪声w(n)的均值和方差，便可由上式生成ARMA(p,q)的数据。

4. 函数说明（arma.h）（关于gauss函数，请查看这里）

#ifndef ARMA_H_
#define ARMA_H_

#include <stdlib.h>
#include "gauss.h"

/*
a —— 双精度实型一组数组，长度(p+1)。ARMA(p,q)模型的自回归系数。
b —— 双精度实型一组数组，长度为(q+1)。ARMA(p,q)模型的平滑平均系数。
P —— 整型变量。ARMA(p,q)模型的自回归阶数。
mean —— 双精度实型变量。产生白噪声所用的正态分布均值μ。
sigma —— 双精度实型变量。产生白噪声所用的正态分布的均方差σ。
x —— 双精度实型一组数组，长度为n。存放ARMA(p,q)模型的数据。
n —— 整型变量。ARMA(p,q)模型数据的长度。
*/
void arma(double a[], double b[], int p, int q, double mean, 
           double sigma, long *seed, double x[], int n);

void arma(double a[], double b[], int p, int q,  
           double mean, double sigma, long *seed, double x[], int n)
{
    int i,k,m;
    double s,*w;
    w=malloc(n*sizeof(double));
    //1.生成高斯序列
    for(k=0; k<n; k++)
    {
        w[k]=gauss(mean,sigma,seed);  //产生白噪声（μ = mean，σ = sigma）
    }
    x[0]=b[0]*w[0];  //计算x[0]
    //2.计算和最大阶数对应数组数据
    for(k=1; k<=p; k++)
    {
        s=0.0;

        //（1）计算AR模型值
        for(i=1; i<=k; i++)
        {
           s+=a[i]*x[k-i];
        }
        //(2)计算MA模型值
        s=b[0]*w[k]-s;
        if(q==0)
        {
           x[k]=s;
           continue;
        }
        m=(k>q)? q:k;
        for(i=1; i<=m; i++)
        {
           s+=b[i]*w[k-i];
        }
        x[k]=s;  
    }
    //3.最大阶数之后的剩余的数据
    for(k=(p+1); k<n; k++)
    {
        s=0.0;
        //（1）计算AR模型值
        for(i=1; i<=p; i++)  
        {
            s+=a[i]*x[k-i];
        }
        //(2)计算MA模型值
        s=b[0]*w[k]-s;
        if(q==0)
        {
            x[k]=s;
            continue;
        }
        for(i=1; i<=q; i++) 
        {
            s+=b[i]*w[k-i];
            x[k]=s;
        }
        free(w);
    }
}

#endif /* ARMA_H_ */

5. 主函数实现（main.c）

#include <stdio.h>
#include "arma.h"
 
int main(void)
{
    int i,n,p,q;
    long seed;
    double mean,sigma,x[200];
    static double a[3]={1.0,1.45,0.6};    //自回归系数
    static double b[3]={1.0,-0.2,-0.1};   //平滑平均系数
    FILE *fp;
    n=200;
    p=2;
    q=2;
    seed=13579l;
    mean=0.0;
    sigma = 0.5;
    arma(a,b,p,q,mean,sigma,&seed,x,n);
    for(i=0;i<32;i+=4)
    {
       printf("    %10.7lf     %10.7lf",x[i],x[i+1]);
       printf("   %10.7lf       %10.7lf",x[i+2],x[i+3]);
       printf("\n");
    }
    fp=fopen("arma.dat","w");
    for(i=0;i<n;i++) {  fprintf(fp,"%d %1f\n",i,x[i]);}
    fclose(fp);
    return 0;
}

6. 执行结果

7. 使用软件QtiPlot，绘制输出波形，如下：

参考：

（1）ARMA模型及其应用     徐静