covar.c

用协方差方法估计AR模型参数

#include<stdlib.h>
void covar(double x[],int n,int p,double a[],double *v,int mode)
{
	int i,j,k,m;
	double cc,sum,*c;
	c=(double *)malloc((p*(p+1)/2)*sizeof(double));
	m=0;
	for(k=1;k<=p;k++)
		for(j=1;j<=k;j++)
		{
			c[m]=0.0;
			for(i=p;i<n;i++)
			{c[m]+=x[i-j]*x[i-k];}
			if(mode==1)
			{
				for(i=0;i<(n-p);i++)
				{c[m]+=x[i+j]*x[i+k];}
			}
			m=m+1;
		}
		for(j=1;j<=p;j++)
		{
			a[j-1]=0.0;
			for(i=p;i<n;i++)
			{a[j-1]-=x[i-j]*x[i];}
			if(mode==1)
			{
				for(i=0;i<(n-p);i++)
				{a[j-1]-=x[i+j]*x[i];}
			}
		}
		cholesky(c,a,p);              //cholesky函数参见cholesky.rar
		for(k=(p-1);k>=0;k--)
		{a[k+1]=a[k];}
		a[0]=1.0;
		sum=0.0;
		for(k=0;k<=p;k++)
		{
			cc=0.0;
			for(i=p;i<n;i++)
			{cc+=x[i]*x[i-k];}
			if(mode==1)
			{
				for(i=0;i<(n-p);i++)
				{cc+=x[i]*x[i+k];}
			}
			if(k==0)
			{sum+=cc;}
			else
			{sum+=cc*a[k];}
		}
		if(mode==1)
		{v[0]=sum/(2*(n-p));}
		else
		{v[0]=sum/(n-p);}
		free(c);
}