模型阶次递推.c

时间序列分析

#include "math.h"
#include "stdio.h"
#include "brmul.c"
#include "brinv.c"

main()
{
	int i,j;
	double z[620],u[620],e[620],zz[600][1],m[600][2],mz[2][600],mm[2][2],t,wg[600][1],p[600][1],pz[1][600], 
	B[2][2],bm[2][600],gg2[2][1],www[2][600],wm[2][2], w1[600][2],w1z[2][600],ww1[2][2],g1[2][1],v1[1][1],
	mw2[2][2],mww2[2][2],wm2[2][2],am2[2][600],ad2[2][1],gg1[2][1],w2[600][4],w2z[4][600],A1[2][2],g2[4][1],
	v2[1][1],www2[2][600],mw3[2][4],mww3[2][4],ww3[4][4],www3[4][600],wm3[4][2],am3[4][600],ad3[4][1],gg3[4][1],
	w3[600][6],w3z[6][600],A2[4][2],g3[6][1],v3[1][1],mw4[2][6],mww4[2][6],ww4[6][6],www4[6][600],wm4[6][2],
	am4[6][600],ad4[6][1],gg4[6][1],w4[600][8],w4z[8][600],A3[6][2],g4[8][1],v4[1][1],mw5[2][8],mww5[2][8],ww5[8][8],
	www5[8][600],wm5[8][2],am5[8][600],ad5[8][1],gg5[8][1],w5[600][10],w5z[10][600],A4[8][2],g5[10][1],v5[1][1],mw6[2][10],
	mww6[2][10],ww6[10][10],www6[10][600],wm6[10][2],am6[10][600],ad6[10][1],gg6[10][1],w6[600][12],w6z[12][600],
	A5[10][2],g6[12][1],v6[1][1];
	FILE *fp1,*fp2,*fp3,*fp4;
	fp1=fopen("mxulie1.txt","r");
	fp2=fopen("baizaosheng1.txt","r");
	fp3=fopen("z.txt","w");
	fp4=fopen("shuchu.txt","w");

	for(i=0;i<620;i++)
	{
	 fscanf(fp1,"%lf",&u[i]);
	   u[i+255]=u[i];
	   u[i+2*255]=u[i];
	 fscanf(fp2,"%lf ",&e[i]);
	}
	z[0]=e[0];
	z[-1]=z[-2]=z[-3]=z[-4]=u[-1]=u[-2]=u[-3]=u[-4]=0;
	for(i=1;i<620;i++)
		z[i]=1.185*z[i-1]-0.814*z[i-2]+0.518*z[i-3]-
			0.349*z[i-4]+0.117*z[i-5]+1.08*u[i-1]-0.745*u[i-2]+
			0.475*u[i-3]-0.253*u[i-4]+0.123*u[i-5]+e[i];
	for(i=0;i<620;i++)
		fprintf(fp3,"%lf  ",z[i]);

	for(i=0;i<600;i++)
	{
		w1[i][0]=-z[i+10];
		w1[i][1]=u[i+10];
		zz[i][0]=z[i+11];
	}
	for(i=0;i<2;i++)
		for(j=0;j<600;j++)
			w1z[i][j]=w1[j][i];
	brmul(w1z,w1,2,600,2,ww1);
	brinv(ww1,2);
	brmul(ww1,w1z,2,2,600,www);
	brmul(www,zz,2,600,1,g1);
	for(i=0;i<2;i++)
	{
		printf("g1[%d]=%lf \n",i,g1[i][0]);
		fprintf(fp4,"g1[%d]=%lf        ",i,g1[i][0]);
	}
	fprintf(fp4,"\n"); 
	brmul(w1,g1,600,2,1,wg);
	for(i=0;i<600;i++)
	{
		p[i][0]=zz[i][0]-wg[i][0];
		pz[0][i]=p[i][0];
	}
	brmul(pz,p,1,600,1,v1);
	printf("v1=%lf\n",v1[0][0]);
	printf("\n");
	fprintf(fp4,"v1=%lf  ",v1[0][0]);
	fprintf(fp4,"\n\n"); 

	for(i=0;i<600;i++)
	{
		m[i][0]=-z[i+9];
		m[i][1]=u[i+9];
	}
	for(i=0;i<2;i++)
		for(j=0;j<600;j++)
		mz[i][j]=m[j][i];

	brmul(mz,m,2,600,2,mm);  
	brmul(mz,w1,2,600,2,mw2);  
	brmul(mw2,ww1,2,2,2,mww2);
	brmul(mww2,w1z,2,2,600,www);
	brmul(www,m,2,600,2,wm);
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
		B[i][j]=mm[i][j]-wm[i][j];
	brinv(B,2);
	brmul(ww1,w1z,2,2,600,www2);
	brmul(www2,m,2,600,2,wm2);
	brmul(wm2,B,2,2,2,A1);
	brmul(A1,mz,2,2,600,am2); 
	brmul(am2,p,2,600,1,ad2);  
	for(i=0;i<2;i++)
		gg1[i][0]=g1[i][0]-ad2[i][0];
	brmul(B,mz,2,2,600,bm); 
	brmul(bm,p,2,600,1,gg2);

	for(i=0;i<2;i++)
		g2[i][0]=gg1[i][0];
	for(i=2;i<4;i++)
		g2[i][0]=gg2[i-2][0];

	for(i=0;i<600;i++)
		for(j=0;j<4;j++)
		{
			if(j<2) 
				w2[i][j]=w1[i][j];
			else 
				w2[i][j]=m[i][j-2];
		}
	for(i=0;i<4;i++)
		for(j=0;j<600;j++)
			w2z[i][j]=w2[j][i];
	brmul(w2,g2,600,4,1,wg);
	for(i=0;i<600;i++)
	{
		p[i][0]=zz[i][0]-wg[i][0];
		pz[0][i]=p[i][0];
	}
	brmul(pz,p,1,600,1,v2);
	for(i=0;i<4;i++)
	{
		printf("g2[%d]=%lf ",i,g2[i][0]);
		printf("\n");
		fprintf(fp4,"g2[%d]=%lf        ",i,g2[i][0]);
	}

	printf("v2=%lf ",v2[0][0]); 
	fprintf(fp4,"v2=%lf ",v2[0][0]);
	fprintf(fp4,"\n"); 
	t=(v1[0][0]-v2[0][0])*(600-2*2)/(v2[0][0]*2);
	printf("t2=%lf\n\n ",t);  
	fprintf(fp4,"t2=%lf ",t);
	fprintf(fp4,"\n\n"); 

	for(i=0;i<600;i++)
	{
		m[i][0]=-z[i+8];
		m[i][1]=u[i+8];
	}
	for(i=0;i<2;i++)
		for(j=0;j<600;j++)
			mz[i][j]=m[j][i];
	brmul(mz,m,2,600,2,mm);  
	brmul(w2z,w2,4,600,4,ww3);  
	brinv(ww3,4);
	brmul(mz,w2,2,600,4,mw3);  
	brmul(mw3,ww3,2,4,4,mww3);
	brmul(mww3,w2z,2,4,600,www);
	brmul(www,m,2,600,2,wm);
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
			B[i][j]=mm[i][j]-wm[i][j];
	brinv(B,2);
	brmul(ww3,w2z,4,4,600,www3);
	brmul(www3,m,4,600,2,wm3);
	brmul(wm3,B,4,2,2,A2);
	brmul(A2,mz,4,2,600,am3);  
	brmul(am3,p,4,600,1,ad3);  
	for(i=0;i<4;i++)
		gg3[i][0]=g2[i][0]-ad3[i][0];
	brmul(B,mz,2,2,600,bm);  
	brmul(bm,p,2,600,1,gg2);
	for(i=0;i<4;i++)
		g3[i][0]=gg3[i][0];
	for(i=4;i<6;i++)
		g3[i][0]=gg2[i-4][0];

	for(i=0;i<600;i++)
		for(j=0;j<6;j++)
		{
			if(j<4) w3[i][j]=w2[i][j];
			else w3[i][j]=m[i][j-4];
		}
	for(i=0;i<6;i++)
		for(j=0;j<600;j++)
			w3z[i][j]=w3[j][i];
	brmul(w3,g3,600,6,1,wg);
	for(i=0;i<600;i++)
	{
		p[i][0]=zz[i][0]-wg[i][0];
		pz[0][i]=p[i][0];
	}
	brmul(pz,p,1,600,1,v3);
	for(i=0;i<6;i++)
	{
		printf("g3[%d]=%lf ",i,g3[i][0]);
		printf("\n");
		fprintf(fp4,"g3[%d]=%lf       ",i,g3[i][0]);
	}
	fprintf(fp4,"\n"); 
	printf("v3=%lf ",v3[0][0]); 
	fprintf(fp4,"v3=%lf",v3[0][0]);
	fprintf(fp4,"\n"); 
	t=(v2[0][0]-v3[0][0])*(600-2*3)/(v3[0][0]*2);
	printf("t3=%lf\n\n ",t); 
	fprintf(fp4,"t3=%lf ",t); 
	fprintf(fp4,"\n\n"); 

	for(i=0;i<600;i++)
	{
		m[i][0]=-z[i+7];
		m[i][1]=u[i+7];
	}
	for(i=0;i<2;i++)
		for(j=0;j<600;j++)
			mz[i][j]=m[j][i];

	brmul(mz,m,2,600,2,mm);  
	brmul(w3z,w3,6,600,6,ww4);  
	brinv(ww4,6);
	brmul(mz,w3,2,600,6,mw4);  
	brmul(mw4,ww4,2,6,6,mww4);
	brmul(mww4,w3z,2,6,600,www);
	brmul(www,m,2,600,2,wm);
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
			B[i][j]=mm[i][j]-wm[i][j];
	brinv(B,2);
	brmul(ww4,w3z,6,6,600,www4);
	brmul(www4,m,6,600,2,wm4);
	brmul(wm4,B,6,2,2,A3);
	brmul(A3,mz,6,2,600,am4); 
	brmul(am4,p,6,600,1,ad4);  
	for(i=0;i<6;i++)
		gg4[i][0]=g3[i][0]-ad4[i][0];
	brmul(B,mz,2,2,600,bm);  
	brmul(bm,p,2,600,1,gg2);

	for(i=0;i<6;i++)
		g4[i][0]=gg4[i][0];
	for(i=6;i<8;i++)
		g4[i][0]=gg2[i-6][0];

	for(i=0;i<600;i++)
		for(j=0;j<8;j++)
		{
			if(j<6) 
				w4[i][j]=w3[i][j];
			else 
				w4[i][j]=m[i][j-6];
	}
	for(i=0;i<8;i++)
		for(j=0;j<600;j++)
			w4z[i][j]=w4[j][i];
	brmul(w4,g4,600,8,1,wg);
	for(i=0;i<600;i++)
	{
		p[i][0]=zz[i][0]-wg[i][0];
		pz[0][i]=p[i][0];
	}
	brmul(pz,p,1,600,1,v4);
	for(i=0;i<8;i++)
	{
		printf("g4[%d]=%lf ",i,g4[i][0]);
		printf("\n");
		fprintf(fp4,"g4[%d]=%lf       ",i,g4[i][0]);
	}
	fprintf(fp4,"\n"); 
	printf("v4=%lf ",v4[0][0]); 
	fprintf(fp4,"v4=%lf ",v4[0][0]); 
	fprintf(fp4,"\n"); 
	t=(v3[0][0]-v4[0][0])*(600-2*4)/(v4[0][0]*2);
	printf("t4=%lf \n\n ",t);
	fprintf(fp4,"t4=%lf ",t);
	fprintf(fp4,"\n\n"); 
 
	for(i=0;i<600;i++)
	{
		m[i][0]=-z[i+6];
		m[i][1]=u[i+6];
	}
	for(i=0;i<2;i++)
		for(j=0;j<600;j++)
			mz[i][j]=m[j][i];

	brmul(mz,m,2,600,2,mm); 
	brmul(w4z,w4,8,600,8,ww5);  
	brinv(ww5,8);
	brmul(mz,w4,2,600,8,mw5); 
	brmul(mw5,ww5,2,8,8,mww5);
	brmul(mww5,w4z,2,8,600,www);
	brmul(www,m,2,600,2,wm);
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
			B[i][j]=mm[i][j]-wm[i][j];
	brinv(B,2);

	brmul(ww5,w4z,8,8,600,www5);
	brmul(www5,m,8,600,2,wm5);
	brmul(wm5,B,8,2,2,A4);
	brmul(A4,mz,8,2,600,am5); 
	brmul(am5,p,8,600,1,ad5); 
	for(i=0;i<8;i++)
		gg5[i][0]=g4[i][0]-ad5[i][0];

	brmul(B,mz,2,2,600,bm);  
	brmul(bm,p,2,600,1,gg2);
	for(i=0;i<8;i++)
		g5[i][0]=gg5[i][0];
	for(i=8;i<10;i++)
		g5[i][0]=gg2[i-8][0];

for(i=0;i<600;i++)
	for(j=0;j<10;j++)
	{
		if(j<8) 
			w5[i][j]=w4[i][j];
		else
			w5[i][j]=m[i][j-8];
	}
	for(i=0;i<10;i++)
		for(j=0;j<600;j++)
		w5z[i][j]=w5[j][i];

	brmul(w5,g5,600,10,1,wg);
	for(i=0;i<600;i++)
	{
		p[i][0]=zz[i][0]-wg[i][0];
		pz[0][i]=p[i][0];
	}
	brmul(pz,p,1,600,1,v5);
	for(i=0;i<10;i++)
	{
		printf("g5[%d]=%lf ",i,g5[i][0]);
		printf("\n");
		fprintf(fp4,"g5[%d]=%lf        ",i,g5[i][0]);
	}
	fprintf(fp4,"\n"); 
	printf("v5=%lf ",v5[0][0]); 
	fprintf(fp4,"v5=%lf ",v5[0][0]);
	fprintf(fp4,"\n"); 
	t=(v4[0][0]-v5[0][0])*(600-2*5)/(v5[0][0]*2);
	printf("t5=%lf \n\n ",t);  
	fprintf(fp4,"t5=%lf  ",t); 
	fprintf(fp4,"\n\n"); 

	for(i=0;i<600;i++)
	{
		m[i][0]=-z[i+5];
	 
		m[i][1]=u[i+5];
	}
	for(i=0;i<2;i++)
		for(j=0;j<600;j++)
			mz[i][j]=m[j][i];

	brmul(mz,m,2,600,2,mm);  
	brmul(w5z,w5,10,600,10,ww6);  
	brinv(ww6,10);
	brmul(mz,w5,2,600,10,mw6);  
	brmul(mw6,ww6,2,10,10,mww6);
	brmul(mww6,w5z,2,10,600,www);
	brmul(www,m,2,600,2,wm);
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
			B[i][j]=mm[i][j]-wm[i][j];
	brinv(B,2);
	brmul(ww6,w5z,10,10,600,www6);
	brmul(www6,m,10,600,2,wm6);
	brmul(wm6,B,10,2,2,A5);
	brmul(A5,mz,10,2,600,am6);  
	brmul(am6,p,10,600,1,ad6); 
	for(i=0;i<10;i++)
		gg6[i][0]=g5[i][0]-ad6[i][0];
	brmul(B,mz,2,2,600,bm);  
	brmul(bm,p,2,600,1,gg2);
	for(i=0;i<10;i++)
		g6[i][0]=gg6[i][0];
	for(i=10;i<12;i++)
		g6[i][0]=gg2[i-10][0];

	for(i=0;i<600;i++)
		for(j=0;j<12;j++)
		{
			if(j<10) 
				w6[i][j]=w5[i][j];
			else 
				w6[i][j]=m[i][j-10];
		}
	for(i=0;i<12;i++)
		for(j=0;j<600;j++)
			w6z[i][j]=w6[j][i];
	brmul(w6,g6,600,12,1,wg);
	for(i=0;i<600;i++)
	{
		p[i][0]=zz[i][0]-wg[i][0];
		pz[0][i]=p[i][0];
	}
	brmul(pz,p,1,600,1,v6);
	for(i=0;i<12;i++)
	{
		printf("g6[%d]=%lf    ",i,g6[i][0]);
		printf("\n");
		fprintf(fp4,"g6[%d]=%lf        ",i,g6[i][0]);
	}

	printf("v6=%lf ",v6[0][0]);
	fprintf(fp4,"v6=%lf  ",v6[0][0]); 
	fprintf(fp4,"\n"); 
	t=(v5[0][0]-v6[0][0])*(600-2*6)/(v6[0][0]*2);
	printf("t6=%lf  \n\n ",t);  
	fprintf(fp4,"t6=%lf  ",t); 
	fprintf(fp4,"\n\n"); 
	fprintf(fp4,"f检验的参数输出结果如下：  \n ");
	for(i=0;i<10;i++)
		fprintf(fp4,"%lf       \n ",g5[i][0]);

	fclose(fp1);
	fclose(fp2);
	fclose(fp3);
	fclose(fp4);

}