lmy.cpp

时间序列分析

 #include <iostream.h>
 #include <fstream.h>
 #include <stdlib.h>
 #include <math.h>
 #define N 900                          //参数递推最后次数
 #define M 300                          //时变参数变化的转折点，M-100为开始递推的起始点
 #define na 1                           //系统阶次
 #define nb 1

 int brmul(double a[],double b[],int m,int n,int k,double c[]);  //矩阵相乘
 int brinv(double a[],int n);                                    //矩阵求逆

 void main()
 {
	double u[N+na+1],e[N+na+1],a[N+na+1],b[N+na+1],y[N+na+1],z[M-100];
	double w1[(M-100)*(na+nb)],w1t[(na+nb)*(M-100)],sita[na+nb],p[(na+nb)*(na+nb)],f[(na+nb)*(M-100)];
	double wn[na+nb],wna[na+nb],sita1[na+nb],p1[(na+nb)*(na+nb)];
	double f1[na+nb],f2[na+nb],f3[1],f4[1],f5[(na+nb)*(na+nb)],f6[(na+nb)*(na+nb)];
	double g1[na+nb],g2[na+nb],g3[1],g4[2],g5[(na+nb)*(na+nb)],g6[(na+nb)*(na+nb)];
	double t1[na+nb]={0.6,0.3},ta[na+nb];
	int i,j,k;
	ifstream fip1("M序列.txt");
	ifstream fip2("Gauss2.txt");
	ofstream fop1("tdata.txt");
	ofstream fop2("mdata.txt");
	for(i=0;i<N+na+1;i++)
	{
		fip1>>u[i];
		fip2>>e[i];
	}
	for(i=0;i<M;i++)
	{
		a[i]=0.8;b[i]=0.5;
		fop1<<a[i]<<"  "<<b[i]<<"  "<<endl;
	}
	for(i=M;i<N+na+1;i++)
	{
		a[i]=0.6;b[i]=0.3;
 		fop1<<a[i]<<"  "<<b[i]<<"  "<<endl;
	}
	y[0]=e[0];
	for(i=1;i<N+na+1;i++)
		y[i]=-a[i]*y[i-1]+b[i]*u[i-1]+e[i];
	//先用M-100组数据一次辩识系统的参数
    for(i=0;i<M-100;i++)
		z[i]=y[na+i+1];
    for(i=0;i<M-100;i++)
	{
		for(j=0;j<na;j++)
			w1[i*(na+nb)+j]=-y[na+i];
		for(j=0;j<nb;j++)
			w1[i*(na+nb)+na+j]=u[na+i];
	}
	for(i=0;i<na+nb;i++)
		for(j=0;j<M-100;j++)
			w1t[i*(M-100)+j]=w1[j*(na+nb)+i];
    brmul(w1t,w1,na+nb,M-100,na+nb,p);
	brinv(p,na+nb);
	brmul(p,w1t,na+nb,na+nb,M-100,f);
	brmul(f,z,na+nb,M-100,1,sita);
	for(i=0;i<na+nb;i++)
	{
		cout<<sita[i]<<"  ";
        fop2<<sita[i]<<"  ";
	}
	cout<<endl;
	fop2<<endl;
	//以上面计算的各参数为初始值进行递推
    for(i=M-100;i<N;i++)
	{
		for(j=0;j<na;j++)
			wn[j]=-y[na+i];
		for(j=0;j<nb;j++)
			wn[na+j]=u[na+i];
		brmul(p,wn,na+nb,na+nb,1,f1);
        brmul(wn,p,1,na+nb,na+nb,f2);
		brmul(f2,wn,1,na+nb,1,f3);
	    brmul(wn,sita,1,na+nb,1,f4);
		brmul(f1,wn,na+nb,1,na+nb,f5);
		brmul(f5,p,na+nb,na+nb,na+nb,f6);
		for(j=0;j<na+nb;j++)
		{
			sita1[j]=sita[j]+f1[j]*(y[na+i+1]-f4[0])/(f3[0]+1);
			for(k=0;k<na+nb;k++)
			{
                p1[j*(na+nb)+k]=p[j*(na+nb)+k]-f6[j*(na+nb)+k]/(f3[0]+1);
			}
		}
        for(j=0;j<na;j++)
			wna[j]=-y[na+i-M+100];
		for(j=0;j<nb;j++)
			wna[na+j]=u[na+i-M+100];
        brmul(p1,wna,na+nb,na+nb,1,g1);
		brmul(wna,p1,1,na+nb,na+nb,g2);
		brmul(g2,wna,1,na+nb,1,g3);
		brmul(wna,sita1,1,na+nb,1,g4);
		brmul(g1,wna,na+nb,1,na+nb,g5);
		brmul(g5,p1,na+nb,na+nb,na+nb,g6);
		for(j=0;j<na+nb;j++)
		{
			sita[j]=sita1[j]+g1[j]*(y[na+i-M+100+1]-g4[0])/(g3[0]-1);
			ta[j]=fabs(sita[j]-t1[j])/t1[j];
				cout<<sita[j]<<"";         //输出递推参数
			fop2<<sita[j]<<"\t";
			for(k=0;k<na+nb;k++)
			{
				p[j*(na+nb)+k]=p1[j*(na+nb)+k]-g6[j*(na+nb)+k]/(g3[0]-1);
			}
		}
        cout<<endl;
		fop2<<endl;
        if(ta[0]<0.01&&ta[1]<0.01)
			break;

	}
	cout<<i-M-100<<endl;

 }

 int brmul(double a[],double b[],int m,int n,int k,double c[])     //矩阵相乘
 { 
	int i,j,l,u;
    for (i=0; i<=m-1; i++)
    for (j=0; j<=k-1; j++)
      { 
		u=i*k+j; c[u]=0.0;
        for (l=0; l<=n-1; l++)
          c[u]=c[u]+a[i*n+l]*b[l*k+j];
      }
    return 0;
  }

 int brinv(double a[],int n)                    //矩阵求逆
  { int *is,*js,i,j,k,l,u,v;
    double d,p;
    is=(int *)malloc(n*sizeof(int));
    js=(int *)malloc(n*sizeof(int));
    for (k=0; k<=n-1; k++)
      { d=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
          { l=i*n+j; p=fabs(a[l]);
            if (p>d) { d=p; is[k]=i; js[k]=j;}
          }
        if (d+1.0==1.0)
          { free(is); free(js); cout<<"err**not inv\n";
            return(0);
          }
        if (is[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=is[k]*n+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (js[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+js[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        l=k*n+k;
        a[l]=1.0/a[l];
        for (j=0; j<=n-1; j++)
          if (j!=k)
            { u=k*n+j; a[u]=a[u]*a[l];}
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
              if (j!=k)
                { u=i*n+j;
                  a[u]=a[u]-a[i*n+k]*a[k*n+j];
                }
        for (i=0; i<=n-1; i++)
          if (i!=k)
            { u=i*n+k; a[u]=-a[u]*a[l];}
      }
    for (k=n-1; k>=0; k--)
      { if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=js[k]*n+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+is[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
      }
    free(is); free(js);
    return(1);
  }