arm.cs

C#编写的AR(M）趋势预测模型函数

using System;
class MathMethod
{
	/// <summary>
	/// 计算得到AR_m模型的系数向量
	/// </summary>
	/// <param name="N">输入时间序列长度</param>
	/// <param name="m">模型阶次</param>
	/// <param name="x">输入时间序列</param>
	/// <returns>系数向量</returns>
	public static double[] ARm(int N,int m,double[] x)
	{
		//求相关系数
		double[] r = new double[m+1];
		for(int i = 0;i < m + 1;++i)
		{
			r[i] = 0.0;
			for(int j = 0;j < N-i;++j)
				r[i] += x[j] * x[j+i];
			r[i] /= N;
		}
		for(int i = 1;i < m + 1;++i)
			r[i] /= r[0];
		//解方程
		double[][] a = new double[m][];
		for(int i = 0;i < m;++i)
		{
			a[i] = new double[m];
			for(int j = 0,k = i;j < i;++j,--k)
				a[i][j] = r[k];
			a[i][i] = 1.0;
		}
		double[] b = new double[m];
		for(int i = 0;i < m;++i)
			b[i] = r[i+1];
		double[] y = new double[m];
		y = Cal(a,b,m);
		return y;	
	}
	/// <summary>
	/// 求解对称方程组的LDLT分解法
	/// </summary>
	/// <param name="a">系数矩阵a，只需给出下三角部分</param>
	/// <param name="n">系数矩阵维数</param>
	/// <param name="n">右端向量</param>
	/// <returns>系数向量</returns>
	public static double[] Cal(double[][] a,double[] b,int n)
	{
		double[] x = new double[n];
		double[] r = new double[n-1];
		//系数矩阵分解
		double sum = 0.0;
		for(int k = 0;k < n;++k)
		{
			for(int p = 0;p < k;++p)
			{
				r[p] = a[p][p] * a[k][p];
			}
			sum = 0.0;
			for(int p = 0;p< k;++p)
			{
				sum += a[k][p] * r[p];
			}
			a[k][k] -= sum;
			if(Math.Abs(a[k][k]) < Double.Epsilon)
				Console.WriteLine("Wrong!");
			for(int i = k+1;i < n;++i)
			{
				sum = 0.0;
				for(int p = 0;p < k;++p)
					sum += a[i][p] * r[p];
				a[i][k] = (a[i][k] - sum)/a[k][k];
					
			}
		}
		//回代
		for(int i = 1;i < n;++i)
		{
			sum = 0.0;
			for(int k = 0;k < i;++k)
				sum += a[i][k] * b[k];
			b[i] -= sum;
		}
		x[n-1] = b[n-1] / a[n-1][n-1];
		for(int i = n-2;i>=0;--i)
		{
			sum = 0.0;
			for(int k = i+1;k < n;++k)
				sum += a[i][i] * a[k][i] * x[k];
			x[i] = (b[i] - sum)/a[i][i];
		}
		return x;
	}
	/// <summary>
	/// 计算预测值
	/// </summary>
	/// <param name="x">原始输入时间序列</param>
	/// <param name="a">AR模型系数</param>
	/// <param name="Num">预测时间序列长度</param>
	/// <returns></returns>
	public static double[] Pre(double[] x,double[] a,int Num)
	{
		double[] y = new double[Num];
		int m = a.Length;
		int N = x.Length;
		for(int i = 0;i < Num;++i)
		{
			y[i] = 0.0;
			for(int j = 0;j < i && j <  m;++j)
				y[i] += a[j] * y[i - 1 - j];
			for(int j = i;j < m;++j)
				y[i] += a[j] * x[N - 1 - j + i];
		}
		return y;			
	}
}