gls.c

时间序列分析

#include <stdio.h>
#include <math.h>
#include "brmul.c"
#include "brinv.c"
#define L 800
#define N 600
void main()
{ int i,j,k;
  double  u[L],e[L],y[L],y1[N],z[N][1];
  double b[3]={-0.5,0.5,1.0},e1[L],v[N];
  double w[N][3],wt[3][N],wtw[3][3],p[3][N];
  double Q[N][1],Qt[1][N];
  double qtq,s=0;
  double ys[L],zs[N][1],us[N],ws[N][3],wst[3][N];
  double wstws[3][3],p1[3][N];

  FILE *fp1,*fp2,*fp3;
  fp1=fopen("m.txt","r");
  fp2=fopen("whitenoise.txt","r");
  fp3=fopen("reslult.txt","w");
  for(i=0;i<L;i++)
  { 
	  fscanf(fp1,"%lf",&u[i]);//读入输入数据和白噪声
	  fscanf(fp2,"%lf",&e[i]);
  }
  fclose(fp1);
  fclose(fp2);
 
  y[0]=e[0];
  y[1]=e[1];
	for(i=2;i<L;i++)//令s为0，求出输出y 
	  y[i]=-b[0]*y[i-1]-b[1]*y[i-2]+b[2]*u[i]+e[i];
    for(i=0;i<N;i++)
	  z[i][0]=y[i+2];

   //------------------设S为0，估计A
  for(i=0;i<N;i++)
  {   for(j=0;j<2;j++)
  {   w[i][j]=-y[i+1-j];
      w[i][2]=u[i+2];
  }
      for(j=0;j<3;j++)
		  wt[j][i]=w[i][j];
  }
  brmul(wt,w,3,N,3,wtw);
  brinv(wtw,3);
  brmul(wtw,wt,3,3,N,p);
  brmul(p,z,3,N,1,b);
  printf("当s为零时估计的A:\n");
  for(i=0;i<3;i++)
  {  printf("%lf\t",b[i]);
     fprintf(fp3,"b[%d]=%lf\t",i,b[i]);
	 }
  fprintf(fp3,"\n");
  printf("\n");
  printf("\n");
  
  y1[0]=0;
  y1[1]=0;
  y1[2]=0;
  for(i=3;i<L;i++)
	  y1[i]=-0.35*y1[i-1]-0.075*y1[i-2]-0.425*y1[i-3]+u[i]+0.85*u[i-1]+e[i];
	for(k=0;k<7;k++)//进行循环
 {   //--------------------估计参数s
   e1[0]=0; 
   e1[1]=0;
  for(i=2;i<L;i++)
       e1[i]=y1[i]+b[0]*y1[i-1]+b[1]*y1[i-2]-b[2]*u[i];
  for(i=0;i<N;i++)
  {	  v[i]=e1[i+2];
      Q[i][0]=-e1[i+1];
	  Qt[0][i]=Q[i][0];
  }
  qtq=0;
  for(i=0;i<N;i++)
	  qtq=qtq+Q[i][0]*Qt[0][i];
    s=0;
  for(i=0;i<N;i++)
	  s=s+Qt[0][i]*v[i];
  s=s/qtq;
  printf("第%d次:\t",k);
  printf("s=%lf\t",s);
  fprintf(fp3,"第%d次：s=%lf\t\t",k,s);
   //-------------改进的估计
  ys[0]=0;
  us[0]=0;
  for(i=1;i<L;i++)
  {  ys[i]=y1[i]+s*y1[i-1];
     us[i]=u[i]+s*u[i-1];
  }
  for(i=0;i<N;i++)
	  zs[i][0]=ys[i+3];
  for(i=0;i<N;i++)
  {	  for(j=0;j<2;j++)
	    ws[i][j]=-ys[i+2-j];
	     ws[i][2]=us[i+3];
	  for(j=0;j<3;j++)
	     wst[j][i]=ws[i][j];
	}
  brmul(wst,ws,3,N,3,wstws);
  brinv(wstws,3);
  brmul(wstws,wst,3,3,N,p1);
  brmul(p1,zs,3,N,1,b);
  for(i=0;i<3;i++)
  {   printf("%lf\t",b[i]);
      fprintf(fp3,"b[%d]=%lf\t",i,b[i]);
  }
  printf("\n");
  fprintf(fp3,"\n");
  }
  fclose(fp3);
}