表 1：2005—2008 年巢湖市年负荷量 年份 年最大负荷 (万千瓦) 2005 2006 2007 2008 2009 28.8956 34.7893 37.8492 42.1131 45.9033 巢湖市年负荷量保持了一定规律性的同时，其不确定性也带来了 负荷变化的随机性。灰色系统 GM(1，1)预测指对获取信息量较少、 数据资料年限较短、波动不大且逐年呈上升趋势的问题进行预测的有 效方法。但由于 MG(1，1)模型的解为指数型曲线，一般呈单调递增， 进行超长期预测，往往偏高，因而对波动性大的数据拟合较差，预测 精度较低。 马尔柯夫预测是根据状态之间的转移概率来预测系统未来变化 的，它反映了各状态之间转移的内在规律和各种随机因素的影响程 度。因此，马尔柯夫预测随机波动性较大的问题较为适宜，正好可以 弥补 GM（1, 1)模型的局限性，而 GM (1, 1)又可揭示预测问题时 间序列发展变化的总趋势，用来弥补马尔柯夫预测的局限性，达到取 长补短。两者结合形成的灰色 GM (1, 1)马尔柯夫预测方法，可大大 提高预测的精度。 就负荷预测变化来看，既存在时间序列变化的规律性，也存在众 多因素干扰导致的随机性。利用灰色 GM(1，l)—马尔柯夫预测模型 对房价进行预测，能有效契合房价的变化趋势。 首先建立灰色 GM(1，l)。 设 n} 0X （ ）（k）|k=1,2, ， 为非负原始时间序列数据，进行一次累加 

生成得: X (1) ( ) k  k  i 1  (0) x ( ), i k  1,2, n  确定数据矩阵 B，Q 1 1 1           (1) (1)  x (1) (2)] x (1) (2)  x (1) (3)]  ( n 1)   x (1) ( ) n (1)  1 [ x 2 1 [  2 B             x 1 [ 2 (0) (2) (3) (0)  (0) ( ) n        Q         x x x 计算参数（最小二乘法） A  a     u    ( T B B 1 ) B  T Q X (1) ( k 1)   ( X (0) (1)  u a  ak ) e  u a 确定 GM(1，l)模型 令 ( )Y k X （0）（k+1)则 X (0) ( k 1)   X (1) ( k 1 X k   （ ）） （ ） （k） Y  0 式中 a 为发展系数，其大小反映了数据序列 0X（ ）的增长速度;u 为内生变量;Y(k) 为预测值。 接着以预测值 Y(k)为基准，在其上下两侧作 m 条与之平行的曲线 Yi(k)二 Y(k)+Ai (Ai 为常数，m、Ai 视具体情况而定)，每相邻两条曲线之间的区域为一个状态， 0X（ ） 将芜湖市市房地产现状分析及未来走势预测符合马氏链的非平稳随机序列{ (k+1)}划分 m 个状态，一记为 1E 2E ，…， mE 。 最后确定马尔柯夫状态转移概率矩阵。数据序列由状态 Ei 经过 w 步转移状 态 Ej 的概率称作 w 步状态转移概率，即 

) P w ij (  ( ) m w ij m ( ) i i , j  1,2,3  m 式中 ( ijm w 为状态 Ei 经过 w 步转移到状态 Ej 的原始数据样本数，Mi 为状 ) 态 Ei ) ) w 的原始数据样本数。则状态转移概率矩阵为: p p （w） 12 1 p p 22 2 m    p p  ( p 11  ( p  21      p     ( 1 m P m m w w w ) w ) w ) 2 ( ( ) ( ( w ) ( w ) ( w ) mn        在利用状态转移概率矩阵分析实际问题时，一般只需考虑一步转移概率矩阵 (1)p 假定预测对象处于 Ei 状态，考察转移概率矩阵 (1)p 中第 i 行，若 max( p ( ij ) (1) ) 则系统未来时刻最有可能从状态 Ei 转向状态 E，。 根据己确定的系统未来时刻状态状态，预测结果最可能为: p （ ）， 1 ij ' ( Y k  ） Y( ) k  1 2 ( A i 上  下 i A ） 建立灰色 GM(1，1)模型 计算 X(1)的紧邻均值生成序列得到计算结果如下： 1 1  1            B               1 2 1 2 1 2          Y N          X X X (1) (1) X X (1)  X (1) ( 2 ) (3) ( 2 )  X (1)        (1) X ( n 1)   X (1) ( n )    ( 0 ) ( 2 ) ( 0 ) ( 3 )  ( 0 ) ( n )         。 对参数列 ˆ a  T , a b    进行最小二乘估计，得 ˆ a  ( T T B B B Y ) 1  。 

即 a=-0.0699，u=2404.1。 (1) dX dt (1)  aX   将 a,u 的值代入 ，确立 GM(1,1)模型的微分方程为： (1) dX dt  0.1033 X (1)  34314.7 则可得时间响应式为： 48671 k 1   ，  ,2,3,4,5 (1) ˆ X ( ) k  (0)    X ( 1) a k    (1) / u a e    / u a  48671e 0.0681 k 1   （ ） 三、灰色程序 function gm1(x) X0=[11.8256 23.7233 32.2392 12.1231 25.1233]; X1=cumsum(X0); N=size(X0); n=2:5; s=n-1; Z1(n-1)=(X1(n)+X1(n-1))/2 ; % Z1=Z1(s); b=[1 1 1 1 ]; B=cat(2,-Z1',b') ; Y1(n-1)=X0(n)'; Y1=Y1(s); Y=Y1'; % a1=inv(B'*B)*B'*Y ; a=a1(1); b=a1(2); k=1:5; X1y(k+1)=(X0(1)-b/a)*exp(-a*k)+b/a; X1y(1)=X1(1) ; % X0y(k+1)=X1y(k+1)-X1y(k); X0y(1)=X1y(1) ; % plot(2005:2010,X0y,'r'); hold on; grid on; plot(2005:2009,X0); hold on; plot(2005:2009,X0,'*'); pause 

a b k=1:18; m(k)=(X0(1)-b/a)*exp(-a*(k-1))+b/a; m; for k=1:18 if k>1 y(k)=m(k)-m(k-1); y(k)=m(k); else end end y