ARMA.rar_ARMA_ARMA 参数 C++_ARMA 模型_ARMA模型_数据模型

#include<math.h> #include<iostream.h> #include<stdio.h> #include<stdlib.h> #include <time.h> #define N 1000 double e[N]; double x[N]; double f0[N][N]={0},Pr[N]={0}; void random(long n,double e[],double avr,double var);//产生高斯白噪声 //产生高斯白噪声 void random(long n,double e[],double avr,double var) { long i; double a,b,c; srand((unsigned)time(NULL)); for(i=0;i<n;i++) { c=rand(); if((c-0.000001)<0) { i--; continue; } a=sqrt(-2.0*log(c/32767.0)); b=2*3.14159265*rand()/32767.0; e[i]=var*a*cos(b)+avr; } } void main () { double Rx[N]={0},r=0; double rx[N]={0}; long i=0,j=0; random(N,e,0,1); //产生N个均值为0，方差为1的高斯白噪声 //得出X[i]的值 x[0]=0; x[-1]=0; //计算X[i]序列函数值 for(i=1;i<=N;i++) x[i]=e[i]-1.5*e[i-1]+0.3*x[i-1]; printf("输出Xt函数值：\n"); for(i=0;i<100;i++) { for(j=0;j<10;j++) { printf("%8.5f ",x[10*i+j]); } printf("\n"); } //计算自相关函数值 for(i=1;i<=N;i++) r=x[i]*x[i]+r; r=r/N; for(i=0;i<N;i++) { for(j=1;j<(N-i);j++) rx[i]+=x[j]*x[j+i]; Rx[i]=rx[i]/(N-i)/r; } printf("输出N/10个自相关函数值：\n"); for(i=0;i<10;i++) { for(j=0;j<10;j++) { printf("%8.5f ",Rx[10*i+j]); } printf("\n"); } //求偏相关函数值 long k=0,p=0,q=0; double sum4=0,sum5=0; f0[1][1]=Rx[0]; Pr[1]=f0[1][1]; for(k=1;k<=N;k++) { sum4=0; sum5=0; for(p=1;p<=k;p++) { sum4=sum4+f0[k][p]*Rx[k+1-p-1]; sum5=sum5+f0[k][p]*Rx[p-1]; } f0[k+1][k+1]=(Rx[k+1-1]-sum4)/(1-sum5); for(p=1;p<=k;p++) f0[k+1][p]=f0[k][p]-f0[k+1][k+1]*f0[k][k+1-p]; } long t=0; for(q=1;q<N;q++) { t++; Pr[q]=f0[t][t]; //最后偏相关函数点存在Pr中 } printf("输出N/10个偏相关函数值：\n"); for(i=0;i<10;i++) { for(j=0;j<10;j++) { printf("%8.5f ",Pr[10*i+j]); } printf("\n"); } //判断符合条件的个数 double B=sqrt(1.0/1000.0); int sum=0,flag1=0,flag2=0; printf("比较数根号N分之一大小为：%32.31f\n",B); for(i=0;i<100;i++)//判断是不是MA模型 if(fabs(Rx[i])<=B) sum++; printf("符合条件的ρk个数为：%d。",sum); if(sum>68.5||sum<67.5) { flag1=1; printf("说明此模型自相关函数不带有截尾性，即不是MA模型！\n"); } else printf("说明此模型自相关函数带有截尾性，即是MA模型！\n"); sum=0;//判断是不是AR模型 for(i=0;i<100;i++) if(fabs(Pr[i])<=2*B) sum++; printf("符合条件的ψkk个数为：%d。",sum); if(sum>95||sum<96) { flag2=1; printf("说明此模型偏相关函数不带有截尾性，即不是AR模型！\n"); } else printf("说明此模型偏相关函数带有截尾性，即是AR模型！\n"); if(flag1&&flag2) printf("说明此模型是ARMA模型！\n"); cout<<"由以上结果可以看出自相关函数和偏相关函数都具由拖尾性，所以可以判断出其为ARMA模型"<<endl; cout<<"对于ARMA模型阶数比较难确定。一般采用由低阶到高阶逐个试探，直到经检验认为模型合适为止。"<<endl; //参数估计 矩法估计参数 //对于较低阶模型，可以采用解方程的方法，直接求出参数估计的表达式 //a1=Rx[1],b1=2*Rx[1]/(1+(1-4*Rx[1]*Rx[1])^0.5), //则，拟合模型为X[i]+a1*X[i-1]=e[i]+b1*e[i-1], double a1=0,b1=0; a1=Rx[1]; b1=2*Rx[1]/(1+sqrt(1-4*Rx[1]*Rx[1])); cout<<"ARMA(1,1)模型为"<<endl; cout<<"X[i]"<<a1<<"*X[i-1]=e[i]"<<b1<<"*e[i-1]"<<endl; //模型检验 //其基本思想是，如果模型正确，则模型的估计值与实际值所产生的残差序列应是随机干扰产生的误差，即其就为白噪声序列，否则，模型不正确。 //检验白噪声，那么只需要看x(t)与x(t+s)的协方差即可：当s=0时等于sigma(x)^2，s不为0时为0，这里的协方差用样本协方差计算； double z1=0,z2=0,m=0; e[-1]=0; for(t=0;t<N;t++) { e[t]=x[t]+a1*x[t-1]+b1*e[t-1]; e[t-1]=x[t-1]+a1*x[t-2]+b1*e[t-2]; z1=e[t]*e[t-1]+z1; z2=e[t]*e[t]+z2; } m=z1/z2; //若m接近于零，则e[t]可作白噪声。e[t]的自相关函数接近等于零。模型适用。 cout<<"m="<<m<<endl; cout<<"e[t]可作白噪声。e[t]的自相关函数接近等于零。所以模型适用。"<<endl; }