基于傅立叶级数牛顿外插法的经济参数预测.doc

2006.4 中国地质大学本科生毕业论文 29第一章 引言1.1 数学趋势预测的目的及意义我国是一个经济持续高速发展的国家，其经济参数的预测对于计划和决策有非常重要的意义。就一个企业而言，管理的好与坏关键在于经营，经营关键在于决策，决策的果断与准确在于预测。我国现已是世界经济发展的火车头，近十几年的经济增长率以9%在高速发展。国际国内的市场瞬息万变，企业要克服产-供-销的盲目性，把生产搞好，提高效率必须开展经济预测。及时掌握市场发展方向，发展趋势，正确了解供求情况，增强科学的预见性，这样才能提高我国企业在国际市场的竞争力，进而才能提高本国的经济效益。目前来说，我国的经济预测居国际领先地位。但对与具有季节性波动变化的经济参数预测很难，主要是季节性的波动与逐年的持续增长在模式上难以分析和描述。本文的方法用傅立叶级数描述波动信息特征，用幂函数描述逐年持续增长趋势，并据此作数学趋势的外推，用以预测下一年的具有波动性特征的参数,其优点承袭原始数据变化趋势，用数学模型来描述未来发展变化，最终达到空间的立体预测效果。本文属于计量经济学中的时间序列分析的内容的方法，现在的时间序列分析方法以曲线拟合和外插法为主。由于课题的方法是纯粹的数值趋势推测，具有独创性，但未做任何机制机理的分析，因此仅作为决策参考，并不能做直接依据使用。课题的重点是数学模型的建立，根据经验数据的波动性和持续增长性运用模型进行描述，并得到合理的结果。此外，经济现象变幻莫测，为了解决外界对经济参数的影响，课题难点在于历史数据的波动处理时，权重系数的设定，应根据不同情况采用几组不同权重系数，对傅立叶系数进行处理,减弱经济参数受不规则因素的影响，达到合理的趋势描述结果。第二章 时间序列波动性特征和持续增长趋势的述2.1 波动特征及傅立叶级数1 波动特征：设时间序列波动性特征即季节变动，指自然条件、生活条件及人们生活习惯影响，社会经济现象在一年内某一个特定时期或就以一年为周期性变化。其特征表现比较稳定且可以预见。如：水果供应的淡季和旺季、某家用电器销量的淡季和旺季、农产品的运输问题。此外，也有人们生活习惯造成的影响。如：春节客运量的增加。那么，波动性的社会经济现象可以用数学模型傅立叶级数来形象的描述它。2 傅立叶级数：我们知道在自然现象或是生产实践中，一些实际问题变化趋势呈周期性。如：某家电企业一年12个月份产量呈周期性波动变化趋势，其产量用函数表示，周期为分别为1月，2月，12月关于的变量。设为周期函数，且周期为，且能展开成傅立叶级数，则： （公式2-1）称为函数的傅立叶级数，称为傅立叶系数。根据三角函数系的正交性经过转换运算后，可得到傅立叶系数：（2-2）实际上周期性函数在物理上往往代表振动与波动的现象中的变化规律，而振动与波动现象在物理上可以看成简谐振动或波的迭加，而在数学上讲就是正弦函数和余弦函数迭加的和，公式的每一项称为一个谐量，如公式（2-3）所示。 （2-3）这种公式的方法称为实用谐量分析，本文将用到此数学表达式来描述经济现象中的波动特征。下面就是的求解问题，此参数可形象反映周期性的波动趋势。我们知道三角函数周期为，用数学来表达就是在周期 中取个数，并求出函数值。换句话来说就是测量数据时总是将周期进行6、8、12等份分来处理。因本文将描述某一年的12个月的周期变化值，因此将周期等分为12。此外等分的份数过小或过大都会产生锯齿波的干扰，从会降低反映周期波动的变动趋势，因此一般取最大值为24,最小不能低于6。具体将在下面章节使用描述，也就是第三章的不等间距点的离散傅立叶变换及拟合系数求解。2.2 时间序列的持续增长趋势与幂函数1 持续增长趋势：设时间序列呈持续增长趋势即指由于某种原因的影响，社会的经济现象在相当长的时间内，逐渐增加或逐渐减少的趋势发展。从另一个方面来讲时间序列在较长的时间内按某中规律稳步增长或下降，或是保持在某一固定趋势，它的变化特征与波动变化恰恰相反。例如：石油的消耗量在某一时间段上升或下降。人口的出生率呈上升趋势，而死亡率呈下降趋势。但另外情形，在长期趋势若干年以后，由于某种外在原因而发生变化。如：新产品的的投入，如果产品的质量好、服务好，它的销量会上升。但因时常的竞争激烈，社会消费的普及率的提高，革新时代的到来，产量就会由上升趋势转为下降趋势。那么，这种时间序列的变化趋势可以用数学模型-幂级数部分的和来形象描述此逐年增长的趋势特征。下面先来介绍常见而应用广泛的一种具有如下形式的函数项级数。2 幂函数：与常数项级数一样，我们把（2-4）称为幂级数的部分和。如果这部分和当n时对区间中的每一点都收敛，那么称级数在区间收敛。此时的极限是定义在区间中的函数，记作:. 这个函数称为级数的和函数，简称和，记作：（2-5）我用其部分和的形式来描述逐年的持续增长趋势，并根据此作数学趋势的外推，即外插法，用以预测下一年具有波动性的参数。3 幂函数求解及表示方法：下面先谈一下幂函数系数的求解问题。设幂函数的表达式为，其数学表达式为：（2-6）Lagerange插值的基函数1（2-7）在算法实现时，先做的幂函数的表达式，例如：当n=2时（2-8）上面的可写成（2-9）在计算机中，幂函数可以由数组来存放各项系数，且下标为对应的指数。即在一个数组中，将按序存储，下标为第项、第项、第2项、第1项、第0项的系数。假设定义的数组为A，在此数组中的存储形式如下表所示：标注3210系数表2-1在这里需要注明：数组一般为浮点数，当计算精度要求比较高的时候用双精度数表示。幂函数的这种表示方法，与代数式的写法完全一致，比较直观，并且容易实现函数的运算，对于具体数值的代入求解也是方便计算的。我们知道插值的目的是求解数学表达式中的系数，本文则利用最小二乘法意义下曲线拟合的方法对最优系数进行求解，在下面的章节，分别介绍当的幂函数系数的求解问题，分别用一次，而二次，三次函数对其持续趋势进行描述，并利用外插法做趋势外推，对未来数据预测。第三章 数学模型的建立3.1 不等距节点的离散傅立叶变换在前面的章节介绍的傅立叶级数，可以形象的描述经济现象的周期性变化，那么会有几种情况会影响方程解出傅立叶系数，其中一种会适合方程求解。即当周期节点为不等距时，如果数据点的个数与傅立叶级数的系数的个数相等，比较适合用方程求出各项傅立叶系数和来。当 时，方程已无法解出系数的值，可以减少傅立叶的个数使之符合，但太小时，波动性的描述将很粗略，无法达到预测效果，无实际意义。当 左右时，一般如果不是可以构造的数据，来源于实际的机关机、物理、化学等方面的数据，其傅立叶的高次项都急剧减小，比、等前几项要小二三个数量级，对数据值的影响不大，所以，一般取最大值为24就足可以描述所要的经济数据了。当时，方程的个数超过未知数的个数，一般情况下，多余的方程为矛盾方程，应该利用最小二乘法意义下的最有解。当时设傅立叶级数如公式（2-1）方程组通式为： （3-1）转换成公式（2-2）形式为：其中为式中要求出的系数，即方程组中的变量。假设已经解出，则代回原方程组中，定义误差为：（）（3-2）误差为的函数。到这里要考虑一下，当时间序列呈波动趋势时，用数学模型：来进行描述。但由于时间序列上的点不可能全在此条波动曲线上，从而会使预测值与实际值存在一定误差。要用一种方法来使得误差最小，这就是下一节要介绍的最下平方法.3.2 最小二乘法意义下的曲线拟合拟合求解当经济现象呈波动趋势的时，要检验这种描述它趋势的数学模型方程是否使误差最小,就是我们下面要谈到的最下二乘法的拟合问题，目的来寻求方程组系数的的最优解,系数的数值反映了经济现象的波动趋势,决定了波形的情况。设误差的平方和为的函数。即满足：(3-3)若使误差的平方和最小，则反映拟合趋势为最佳。一般来说就是要求的最小二乘解。由为微积分的知识可得这一求解可归为极值问题，即满足：代入到（3-3），得到：即满足通式：若将（3-1）式写成矩阵的形式的话，则上式为：（3-4）相当于在原来方程的增广矩阵上左乘了一个系数矩阵的转置矩阵，方阵此时变为个未知数，个方程，这样的结果由矩阵知识得，一般有唯一的解值，这就是最小二乘法意义下的方程组拟合求解，也可称为正规方程求解，这样便解决了当数据个数多于傅立叶个数时，无法求解系数的问题，其解法利用高斯主元素法便于求解，得到结果就是傅立叶级数各项的系数，也就是公式（2-3）中的，将系数值代入到此公式中，解以为周期关的函数，即：（3-5）其函数值就是具有波动性特征的结果,其数值特点反映了波动趋势效果的变化特征。方阵的构成及转换方阵伪代码如下：M1=X(M) & 周期值M2=M1IF M124 M2=12ENDIFm2=8 确定 COS 的级次，从 0 开始 确定 SIN 的级次，从 1 开始 DIME A1(M1,M2+1),A(M2,M2+1),FRI(N+1,M2)A=0FOR L=1 TO N FOR I=1 TO M1 0级次 COS项系数 FOR J=1 TO MC J级次 COS项系数 NEXT FOR J=1 TO MS J级次 SIN项系数 NEXT 每个月份的数值，作为常数项 NEXT IF M2=M1 FOR I=1 TO M2 FOR J=1 TO M2+1 A(I,J)=A1(I,J) NEXT NEXT 生成正则矩阵（当行大于列时） A=0 FOR I=1 TO M2 FOR J=1 TO M2+1 FOR K=1 TO M1 经过转换后形成正则矩阵3.3 傅立叶系数的加权处理我们知道傅立叶级数具有周期的波动性，很容易受到外界因素的干扰，这样的话会使预测的结果产生一定的误差，如：在某一年当中，由于自然灾害、社会的原因会影响某几个月份的产值或产量波动的趋势，从而影响到一年的整体趋势结果。为了解决由外界因素干扰的难题，可利用数学模式上的一些变化来消除影响。方法就是傅立叶级数系数的加权处理。首先，我根据实用谐量分析的内容，求出每一年的平均值，求法如下：（3-6）分别计算每个系数占每一年平均值的比重：即：根据比重值，来设定权重值，原则是每年较平稳变化的权值应设相对小一点，而相反有意外事件的时候，如：印度洋海啸、政局混乱造成的影响，其权值就应相对设大一点，这样的处理就解决了预测参数受外界的影响，达到预测的准确性。但，权重系数的设定有一定的难度，需要综合多方面因素进行估算，这样一来虽然解决了外界因素的干扰，可能会产生到预测值和实际值的偏差，为了解决这一难点，可用到另一种方法将每一项系数归一化处理即：这样就会使的波动性曲线既平稳而有较准确描述出经济现象的波动特性，但它仍然未能完全解决每年经济参数受人为及外界因素的影响。本文的难点就是权值的设定，根据每年的实际情况，给出各年的权值，使之乘以每年的傅立叶系数，即：系数，然后在用经过处理后的系数除以权值之和，得到的系数就是理想中的数值，它结合了数学原理和实际情况，便于使大自然中的经济现象客观地展现在数学模型当中，从原理上解决了外界因素的影响。如图（3.1）所示，蓝色波动曲线为预测趋势结果，不仅继承了历史发展的结果，而且将外界因素也客观的考虑在内，达到了合理的描述结果。图3.1运算过程伪代码如下：FOR L=1 TO N FOR J=2 TO M2 将生成的傅立叶系数除各年的平均值 NEXT JNEXT 权值设为1DO WHILE .T. 加一控件，选择权重系数 K=0 FOR I=1 TO N K=K+P(I) NEXT FOR J=2 TO M2 FRI(N+1,J)=0 FOR I=1 TO N 各年份的傅立叶系数乘权重系数 NEXT 每年的系数除以权重系数之和，生成具有比例权限的傅立叶系数3.4 外插值计算及拟合求解在前面所叙述到Lagrange插值可将一组离散数据转化成公式（2-9）的形式，即三次多项式函数，但是有一个明显的特点，就是它不具有承袭性，即每增加一个节点的时候，不仅要增加求和的项数，而且还得将以前的各项也必须重新计算，这样一来造成了计算量的浪费。为了解决这一缺点，课题改用具有承袭性的Newton插值来进行插值运算。利用这种插值方法既使预测的未来结果具有持续增长趋势，而且也减少了计算量。3.4.1 Newton基函数首先，做出几项多项式的N（x）（3-7）使满足由公式（3-7）得出个多项式，即：（3-8）为了便于记忆，将公式（3-4）简化成：由个不同结点可以唯一确定一个次多项式，称为插值以为节点的基函数，即：从上面运算过程取得基函数，便于求解关于个年份的节点的值，转换成本课题来说就是利用一次，二次，三次函数来描述逐年持续增长趋势，以此来做数学趋势的外推，来预测年的值，即牛顿外插法，如图（3-2）多所示。接下来问题是系数求解，常用的方法是利用差商。但本课题用最小二乘法意义下的多项式曲线拟合来求解系数。图（3-2）3.4.2 多项式拟合求解前面的第一节已经介绍过用最小二乘法求解傅立叶系数，下面继续用此方法求解多项式系数的最优解。给定年的一组数据,求。在前面的第一节已经涉及到了最小二乘法意义下的拟合求解问题，在这里将不加叙述，其求解也将归结为求多元函数的极值问题。下面，将构造多元函数的关于系数的正则方程组。则此时的拟合的系数方程为上式方程可转化为下面形式写成矩阵形式为（3-9）前面已经介绍过正则方程，其解存在且唯一。通过解公式（3-9）关于的正则方程，得到方程的最优解，由于多项式的次数越高，其预测的结果的偏差随之加大，因此采用一次，二次，三次函数预测其结果,并根据趋势预测未来值。其构造正则矩阵的伪代码如下：FRI(N+1,1)=1 Q=1 S=0 S(1)=N FOR I=2 TO 7 构成正则方程的系数矩阵 NEXT NEXT T=0 FOR J=1 TO N 构成正则方程的增广部分3.5 线性方程组求解线性方程组求解有很多种方法，在科技、工程、医学、经济等各个领域中，很多问题常常归结为线性方程。一些问题的数学模型虽不直接表现为线性方程组，可其数值解法却需要将该问题“离散化”或“线性化”为线性方程组。如：电学中的网络问题，经济学中的投入产出问题，用最小二乘法求实验数据的曲线拟合问题，工程项目中的三次样条函数的插值问题，用迭代法解非线性方程组的问题，用差分法或者有限元法解微分方程问题等都导致求解线性代数方程组。高斯消去法是一种古老的方法，基于高斯消去法的基本思想而改进、变形得到的住元素消去法、三角分解法仍然是目前计算机上常用的有效方法。本课题将使用高斯住元素消取法对方程组进行求解。3.5.1 高斯主元素消去法1 使用高斯主元素消去法的条件：如，方程能用高斯消去法求解的充分必要条件是的各阶顺序主子式均不为零。若消元过程中允许对增广矩阵进行行交换，则方程组可用消元法求解的充要条件是可逆。2 高斯列主元素算法实现：目的：解线性方程组 ，其中，,。组成：主要包括四个环节：（1） 选主元素；（2）换行；（3）消元；（4）回代。3 存储方式及算法步骤：用二维数组按行存放系数矩阵的元素，用一维数组存放常量的各元素。计算过细后中，采用紧凑存储方式以节省内存。在第步消元时注意到计算，以不必再保存。因此将存入内。同理，将，回代时将存入到。步1 按列选主元素步2 若则输出（为奇异）：停机。步3若，则转步4；否则换行：步4 计算乘子步5 消元计算步6 回代求解以四阶方程组为例列出矩阵行列式如下 增广矩阵为对增广矩阵转换成梯阶阵，由第一列起寻找是否存在绝对值（）大于的情况，若存在（设）则交换T第一行和第三行数据得到矩阵高斯列主元素求解过程程序如下：FUNC LineEqua(C,N) & 线性方程组求解函数，高斯主元素法,N 为未知数个数，解在第 N+1 列返回PRIVATE I,J,K,M,TFOR M=1 TO N T=ABS(C(M,M) K=M FOR I=M+1 TO N & 寻找主元素 IF TABS(C(I,M) T=ABS(C(I,M) K=I ENDIF NEXT IF MK & M列发现绝对值更大的元素，应交换到主元素行 FOR J=M TO N+1 T=C(M,J) C(M,J)=C(K,J) C(K,J)=T NEXT ENDIF IF ABS(C(M,M)1E-7 & 方程无解，返回，N=0 为标志 N=0 C(1,1)=-1E20 RETURN ENDIF FOR I=M+1 TO N & 消元 T=-C(I,M)/C(M,M) FOR J=M+1 TO N+1 C(I,J)=C(I,J)+T*C(M,J) NEXT NEXTNEXTC(N,N+1)=C(N,N+1)/C(N,N) & 得到第 N 个 未知数的值FOR I=N-1 TO 1 STEP -1 & 回代 FOR J=I+1 TO N C(I,N+1)=C(I,N+1)-C(I,J)*C(J,N+1) NEXT C(I,N+1)=C(I,N+1)/C(I,I)NEXT第四章 模型应用如下表（4-1）所示，数据的特点既有每年周期性波动又有逐年增长趋势的特征，分别对其各自特点采用两中数学模型来进行描述，即：傅立叶级数和幂函数。下面通过已知数据来进行模型应用。表4-1月 年 1234567891011122002100102120117105122134139121114109105200311011113213011613314715413312512011520041201201441431191451601711461341301282005132130160157130160175189161146142140.1 波动趋势描述：经过前面的算法，便可求得傅立叶级数的各项系数，代入到公式（3-5）中可得到波形曲线公式，结果如下：（4-1）课题采用最小二乘法意义下的曲线拟合方式，其与原始数据波动趋势逼近，但公式与客观实际情况的近似程度如何，还要靠实践检验。在这里为了削减傅立叶系数对整个波动趋势的影响，可将每个傅立叶系数用数学意义上的处理，在前面已经介绍的处理方法，即将各项傅立叶系数近似视为1（归一化处理），经过处理完成最佳的波动趋势效果。但,外界的变化是不可预测的,因此还要进一步对系数进行处理,也就是前面所说的权重值的手工设定,对其不规则周期波动问题进行解决。根据不同的情况对每一年的傅立叶系数进行权重系数的设定， 这种处理效果属于动态的原理性的设定,因此应结合不同的发展趋势拟定权重系数,进行定量形式化及结构描述。上面的公式（4-1）是按理想状态下描述的结果，未将人为及自然因素考虑在内，前面已经介绍过将系数通过数学意义上的转换成，这样便达到了整体定量形式化的描述，完成系统性的最佳描述效果。2 预测结果：通过前面的算法，将分别求出一次，二次，三次函数的系，代入到公式（2-4）中得到下面公式：代入周期求值，据此持续增长趋势外推，预测年的经济参数，然后将预测值做傅立叶逆变换生成各项傅立叶系数如公式（2-2），在代入到公式（4-1）中，生成各月产值，其预测结果既反映了逐年的持续增长趋势又描述了一年的周期波动特征，又完成了对此经济现象的综合描述并预测出了下一年具有波动性的经济参数。如下表所示针对原始数据趋势外推预测，分别用一次，二次，三次函数预测出2006年不同的趋势结果。表（4-2）月 年 123456789101112200614114117016714617118819917115915314920061431441721691461741912021741611551512006147148177174150178196208179166160155从上表的趋势预测结果可以看出，数学表达式的次数越高，其数值的变化幅度就越大，准确性就会下降，因此不能采用高次数的表达式对持续增长趋势进行描述并预测。第五章 总结显然,相对于复杂的波动性经济现象，仅有定性描述是远远不够的，还需有定量形式化及结构描述。此经济参数预测数学模型，作为人工系数设定的数学原理，无疑是理想的。它不仅能从变化发展上反映现象的整体活动特征，而且能从有机整体联系上,描述事物周期性的关系和结果，且符合经济现象变化的整体特征。虽也是系统思想的一种整体定量形式化描述，同样形式化地说明了经济现象整体发展变易的基本规律，为经济预测领域提供一种认识外界因素影响的整体模式，表达了数学趋势描述客观事物的一些本质特性，虽然此数学模型与其它数学模型相比,则解决了季节性的波动与逐年的持续增长在模式上的难以描述和分析的难点，并据此数学趋势预测了未来数据，利于决策者在未来的几年做出经济战略调整，使企业能处于优势地位，但自然与人为因素的影响是无法预测的，因此，本课题还存在着机制机理方面的不足，仅做参考，不能作为直接的决策依据。 此数学模型立足整体统筹全局，使整体和部分辩证地统一起来；使质变与量变、有限与无限、及结构与功能的考察内在地统一起来；同时使复杂的经济现象在与真实结构相似同构的基础上，形式化地展现于三维空间；体现了宇宙演化必然性与偶然性的统一；为复杂外界变化的运动状态和组织的把握, 提供了静态描述依据；为复杂系统的人工模拟提供了有效的数学原理。致谢整个毕业设计过程，由于自己的在有些知识上掌握不牢固和理解深度的不足，所以真的要感谢我的指导老师对自己的指导和帮助，在做毕业设计的期间，督促指导自己完成毕业设计，特别是在平时朱老师一有时间就询问我们的论文情况，耐心的给我们讲解每个同学在设计当中遇到的问题,更提出了对论文建设性的意见,为我顺利、高效完成毕业设计提供了条件。此外还要特别感谢自己身边的同学，在自己的程序编码和最后调试阶段，给了自己许多灵感与帮助，从中我体会到了团队精神的伟大之处。最后还要感谢一直关心和爱护我的家人，是他们给了自己机会能够完成自己的学业，能够在自己的黄金年龄学到专业的科学文化知识。是家人给了我学习的机会，是家人给了我精神的支柱，让我在困难面前没有被打倒，养成了个人勇于挑战困难的自信与决心。参考文献1 崔国华 计算方法 华中科技大学出版社2 李正龙 经济预测与决策方法 安徽大学出版社3北京大学数学力学系与数学分析与函数论教研室 数学分析（上册） 人民教育出版社4史济民 汤观全 Visual FoxPro 及其应用系统开发 清华大学出版社5南京工学院数学教研组 积分变换 （第三版）高等教育出版社6黄益平 宋立刚 应用数量经济学 上海人民教育出版社7黄铎 陈兰平 王凤 数值分析 科学出版社附论文主程序& 经济参数预测主程序SET TALK OFFSET SCOR OFFSET STAT OFFSET SAFE OFFCLEAR ALLCLOSE ALLCLEARSET DECI TO 12& SET PROC TO PR1& SET PROC TO *PI=3.141592653589PI2=PI*2SELE 1USE 原始数据001M=FCOUNT()-1 &搜索记录N=RECC()-1IF M6 12,10 SAY 数据的列数小于6，无法做傅立叶分析，按任意键退出 WAIT ENDIFIF N24 M2=24ENDIFm2=8MC=INT(M2/2) & 确定 COS 的级次，从 0 开始 MS=M2-MC-1 & 确定 SIN 的级次，从 1 开始 DIME A1(M1,M2+1),A(M2,M2+1),FRI(N+1,M2)A=0FOR L=1 TO N & 建立线性方程组，求解傅立叶系数 FOR I=1 TO M1 A1(I,1)=1 & 0级次 COS项系数 FOR J=1 TO MC A1(I,J+1)=COS(J*X(I)/M1*PI2) & J级次 COS项系数 NEXT FOR J=1 TO MS A1(I,J+MC+1)=SIN(J*X(I)/M1*PI2) & J级次 SIN项系数 NEXT A1(I,M2+1)=C(L,I) & 常数项 NEXT IF M2=M1 FOR I=1 TO M2 FOR J=1 TO M2+1 A(I,J)=A1(I,J) NEXT NEXT ELSE & 生成正则矩阵，当行大于列时 SUSP A=0 FOR I=1 TO M2 FOR J=1 TO M2+1 FOR K=1 TO M1 A(I,J)=A(I,J)+A1(K,I)*A1(K,J) NEXT NEXT & FOR K=1 TO M1 & A(I,M2+1)=A(I,M2+1)+A1(K,I)*C(L,K) & 常数项 & NEXT NEXT ENDIF LineEqua(A,M2)& 调线性方程组求解系数 FOR J=1 TO M2 FRI(L,J)=A(J,M2+1) & 线性方程组的解级即傅立叶系数写入数组 NEXTNEXTFOR L=1 TO N & 归一化处理 FOR J=2 TO M2 FRI(L,J)=FRI(L,J)/FRI(L,1) NEXT JNEXTYEAR1(N+1)=YEAR1(N)+1 & 生成年份标识YEAR1(N+2)=YEAR1(N+1)+.2 & 整数部分表示年份YEAR1(N+3)=YEAR1(N+1)+.3 & 小数部分表示拟合级次YEAR1(N+1)=YEAR1(N+1)+.1suspP=1DO WHILE .T. & 选择权重系数 & 加一控件，生成权重系数 & K=0 FOR I=1 TO N K=K+P(I) NEXT FOR J=2 TO M2 FRI(N+1,J)=0 FOR I=1 TO N FRI(N+1,J)=FRI(N+1,J)+FRI(I,J)*P(I) & 各年份傅立叶系数乘权重系数 NEXT FRI(N+1,J)=FRI(N+1,J)/K NEXT FRI(N+1,1)=1 &平均值a0置为1 Q=1 S=0 S(1)=N FOR I=2 TO 7 & 正则方程系数 FOR J=1 TO N Q(J)=Q(J)*J S(I)=S(I)+Q(J) NEXT NEXT T=0 FOR J=1 TO N & 正则方程常数项 T(1)=T(1)+FRI(J,1) T(2)=T(2)+FRI(J,1)*J T(3)=T(3)+FRI(J,1)*J*J T(4)=T(4)+FRI(J,1)*J*J*J NEXT FOR L=2 TO 4 & 外插值计算，T(5)=线性插值结果, T(6)=抛物线插值结果, T(7)=三次插值结果 FOR I=1 to L & 生成正则方程，将Si,Ti写入数组A矩阵 FOR J=1 TO L A(I,J)=S(I+J-1) NEXT A(I,L+1)=T(I) NEXT LineEqua(A,L) & 用最小二乘法求曲线拟合结果 T(L+3)=0 & 计算外插值,分别用1,2,3次函数预测第 N+1 年的平均值 K=1 FOR J=1 TO L T(5)=T(5)+K*A(J,L+1) K=K*(N+1) NEXT NEXT FOR J=1 TO M2 FRI(N+1,J)=FRI(N+1,J)*T(5)& 傅立叶系数乘预测的平均值 NEXT FOR J=1 TO M1 & 做傅立叶逆变换，生成各月值 C(N+1,J)=FRI(N+1,1) FOR I=1 TO MC & 计算COS 项 C(N+1,J)=C(N+1,J)+FRI(N+1,I+