数学建模提高班专题5——时间序列建模.ppt

2015数学建模提高班差分方程与时间序列模型专题 梦想点燃激情 激情成就未来 柴中林2015 5 9 课件提纲 1差分方程的引例与概念2特殊差分方程的解3平衡点及其稳定性4差分方程组5时间序列与其中的趋势分析6自回归模型7自回归模型识别及参数确定8自回归模型预测及相关说明9建模练习题 1差分方程的引例与概念 例1 某人贷款80万元买了一套房子 期限20年 已知贷款月利率为5 请问他每月要还贷多少 在高数中 我们研究的函数中的变量的取值大都是连续的 在连续区间上取值 如 1 12 但在经济管理和很多实际问题中 变量只能取1 2 3 这样的值 这些变量称为离散型变量 描述离散型变量间关系的模型称为离散型模型 差分方程就是常见的一种离散型模型 微分方程 连续变量间存在函数关系 知道了这个关系 就能够研究变量间的联系与变化规律 然而 这个关系是未知的 但我们可以建立起含自变量 因变量及其导数或微分的等式 这就是微分方程 通过对方程的研究以求得这个函数关系 或通过方程直接揭示变量间的联系就构成了微分方程的主要研究内容 差分方程与微分方程是类似的 只是这里的变量是离散的 差分方程 含自变量 未知函数 因变量 未知函数差分的等式 建立差分方程 求解它的目的就是研究离散变量间的关系 一般的 对有函数关系的两个变量 常用x当自变量 y当因变量 但在差分方程中 因自变量只取整数值 如0 1 2 我们更喜欢用n 或t 表示自变量 这时因变量可用x或y表示 其函数关系是x x n 但我们更常用xn表示 当然 这个关系是不知道的 但我们常能得到的是如下的式子F n xn xn 1 xn k 0 1 或G n xn xn 1 xn k 0 2 或H n xn xn kxn 0 3 这种式子就是差分方程 有时 1 也写成如下的形式xn f n xn 1 xn k 4 因此 差分方程也称为递推关系 考虑例1 用n表示月份 n 0表示贷款月份 xn表示第n月还贷后还欠的钱 r a分别表示银行月利率和月还钱数 xn表示了账户中欠钱数与月份间的函数关系 未知 但我们容易得到一个式子xn rxn a xn 1即xn 1 1 r xn a 5 此外 还有初始条件x0 80 万元 及x240 0 这就是贷款问题的差分方程模型 变化建模比较微元法 对离散关系xn 其函数值构成序 数 列 xn x1 x2 x3 记 xn xn 1 xn 序列后项减前项构成的序列 称为xn的一阶差分 2xn xn xn 1 xn xn 2 2xn 1 xn称为xn的二阶差分 依次类推 对式 5 xn 1 1 r xn a 也可将它写为 xn 1 r xn 1 a或 xn rxn a 差分方程因此而得名 即同一个关系用不同视角不同符号式子会不同 但可以互化 它们是同一个东西 差分方程中的最高阶差分的阶或因变量的最大下标与最小下标之差称为差分方程的阶 差分方程的解是函数 通常有无穷多个 通解是全部解的集合 体现在独立任意常数上 其个数与方程阶数相同 另外 在实际问题中常会给出一些附加条件 如x0的值 称为初始条件 满足初始条件的具体的解就是特解 差分方程问题的研究内容 1差分方程的建立 离散变量关系的建立 也可将连续问题离散化 2差分方程的求解和分析 差分方程在实际问题中有广泛的应用 差分方程的求解并不比微分方程容易 大部分差分方程是无法求解的 这里介绍最简单同时用处很大的一类特殊差分方程的求解 常系数线性齐次差分方程 其一般形式为xn a1xn 1 akxn k 0 6 其中a1 ak是常数 方程 6 有解 其求解步骤为 步骤1 求解对应的特征方程 k a1 k 1 ak 0 7 步骤2 根据步骤1的解的情况写出 6 的通解 2特殊差分方程的解 情况1 若 是 7 的一个单实根 则 n是 6 的一个特解 若 1 2 k是 7 的k个全部不同的单实根 则 6 的通解为xn C1 1n C2 2n Ck kn C1 C2 Ck是任意常数 情况2 若 是 7 的k重实根 则 n n n nk 1 n都是 6 的特解 情况3 若 i是 7 的单重复根 则 ncosn 与 nsinn 都是 6 的特解 其中 是 的模与幅角主值 情况4 若 i是 7 的k重复根 则 ncosn n ncosn nk 1 ncosn 与 nsinn n nsinn nk 1 nsinn 都是 6 的特解 其中 是 的模与幅角 最后 将各个特解如情况1那样与任意常数相组合就得 6 的通解 常系数线性非齐次差分方程 其一般形式为xn a1xn 1 akxn k b n 0 8 8 的求解方法是先求相应齐次方程的通解 记为xn 再求 8 的一个特解 记为xn 0 方法 根据b n 的特点将xn 0 的形式设出 再用待定系数法确定其中的系数 于是 8 的通解为xn xn xn 0 此外 不同于微分方程 对差分方程 当初始条件给定后 可迭代求得任意xn的 精确 值 从而可以对xn的变化规律进行作图分析 如对方程xn f n xn 1 xn k 若x1 x2 xk给定 就可以根据方程依次算出xk 1 xk 2 xk 3 来 下面求解例1 xn 1 1 r xn a 它是一阶常系数非齐次线性差分方程 先解相应的齐次方程xn 1 1 r xn 特征方程为 1 r 其通解为xn C 1 r n C为任意常数 再求其一个特解 从方程看设xn为常数 记为x 代入得xn 0 a r 于是得方程通解 xn C 1 r n a r 代入初始条件得方程组 解之得 大约是5731元 3平衡点及其稳定性 差分方程虽可用迭代法进行数值计算 但计算总归只能进行有限步 其深层次的性质必须用其它工具进行分析 平衡点就是其中一个 平衡点相当于稳定点或不动点 对方程xn f n xn 1 xn k 来说就是若xn 1 xn k都取某一常数 比如a 那么xn也一定是a 从而xn 1 xn 2 xn 3 也都将取a 平衡点就是所有xn都取相同的值 且能使方程成立的点 于是将xn f n xn 1 xn k 中所有xn都换成x 得方程x f n x x 将其求解 每一个解就是一个平衡点 设a是方程的一个平衡点 xn是方程的任一解 若总有则称a是差分方程的一个稳定的平衡点 为什么 稳定的平衡点在实际问题中有重要的价值 现考虑差分方程xn a1xn 1 akxn k 0 并且其解是如下形式xn C1 1n C2 2n Ck kn 显然0是方程的一个平衡点 不难发现对任意s若有 s 1 则必有这说明0是稳定的平衡点 这也是一般差分方程平衡点稳定性的判别方法 若齐次方程的特征方程的根的绝对值都小于1 平衡点稳定 而若某个 的绝对值大于1 平衡点不稳 当 等于1时 有多种情况且实际意义不大 不做讨论 若特征根是复根 就用其模来判断 例2考虑数学模型书中供需关系的蛛网模型 xk 第k时段商品数量 yk 第k时段商品价格 需求函数yk f xk 供需平衡点为P0 x0 y0 当商品生产者的生产只盯着前一期价格 供应函数为xk 1 g yk 时 在平衡点附近各时段商品数量的差分方程模型为xK 1 xk 1 x0 其齐次方程的特征方程的特征根为 所以 1就稳定 否则就不稳 而当商品生产者的生产同时盯着前面两期的价格 供应函数为xk 1 g yk yk 1 2 时 在平衡点附近各时段商品数量的差分方程模型为2xK 2 xk 1 xk 2 1 x0 其齐次方程的特征方程为2 2 0 特征根为 当 8时 根为实根 必有一根绝对值大于1 当0 8时根为复根 用复数的模来判断 可以得到当0 2时稳定 否则不稳 差分方程组 自变量一个 因变量多个 仅讨论线性 线性差分方程组的一般形式为其中aij和bi i j 1 2 n 都是常数 4差分方程组 令记x t x1 t x2 t xn t T b b1 b2 bn T 则上述方程可记为x t 1 Ax t b 该式类似于前面的一阶常系数线性差分方程 可编程数值计算分析 也可利用线性代数理论 主要是特征值和特征向量 进行分析讨论 若x 向量 是该方程的一个平衡点 即x Ax b 则它稳定的条件是A的所有特征值的绝对值都小于1 若某一个的绝对值大于1 就不稳 5时间序列与其中的趋势分析 时间序列 按时间 有时是长度或温度 顺序排列的随机变量序列 但在应用中又指将某个统计指标在不同时间上的各个数值 按时间先后顺序排列而形成的序列 一般等间隔 时间序列分析 根据观测得到的时间序列数据 其机理未知 通过曲线拟合和参数估计来建立数学模型和理论 希望从中寻找出变量的变化规律 对未来的某些阶段进行预测 时间序列有广泛的应用 设 yt 是时间序列 虽然它暗含了时间变量t 但它仅指采样的时间点 因此 一般的不能认为y是纯t的函数 从而按回归等其他理论去做 因为许多变量都随着时间的变化而变化 所以时间序列中也常常包含因时间变量而产生的趋势变化 另外 在时间序列中 相近的各项间往往有很强的依赖关系 当前的数值对下面的数值有很强的影响 如股市 期货 此外 每个数据还受到无法刻画捕捉的随机因素的影响 通常yt可表为yt f t xt 其中f t 表示随时间变化的确定性趋势 xt则主要由随机因素或其积累而形成 是一个平稳序列 在yt f t xt中 趋势成分f t 起着主导的作用 当它存在时 xt可以认为是随机误差 并予以忽略 故可以用回归方法确定f t 中的参数 得到f t 影响f t 的因素有长期趋势 季节变动 季节性规律作用产生的周期变化 循环变动 周期长短不固定的一种变化 以及不规则的变动等 通常 趋势成分主要讨论长期趋势和季节变动趋势 这里也是 当f t 是由长期趋势决定的 其表达式可能是线性趋势f t a bt 二次曲线趋势f t b0 b1t b2t2或更高阶多项式趋势幂函数曲线趋势f t atb对数曲线趋势f t a blnt双曲线趋势f t a b t 或1 f t a b t指数曲线趋势f t aebt修正指数曲线趋势f t L aebt 或f t L abt a 0 0 b 1 龚泊兹曲线趋势 0 a 0 0 b 1 皮尔曲线趋势f t L 1 ae bt 季节规律应是f t s f t 其中s是季节长度还有季节规律与上述趋势的结合模型 f t 应该是什么样的结构 可作图 根据图形的特点进行选择 也可结合事物本身的特点去考虑 有时 最初可选择多个看起来合理的模型 根据模型特点将其转化为线性模型 利用回归等拟合确定其中的参数 再根据得到的模型根据已有的的数据对模型进行评价 从总体误差和未来趋势的角度 最终确定一个最合适的模型 例3下表是1952 1983年我国的社会商品零售总额 对表中数据 做图如下 由图可知 零售额有明显的非线性增加趋势 故趋势曲线可用二次函数 三次函数或指数函数拟合 不能用插值 因为零售额受随机因素的影响 让曲线一定过给定点不适合 例4下表是1961 1981年我国搪瓷脸盆销售数据 对表中数据 做图如下 由图和实际问题可知 销量在中间有显著的增长 之后增长变的缓慢并停顿下来 故趋势曲线可用龚伯兹曲线或皮尔曲线函数 有上限 S型曲线 确定模型中参数对例3 要确定模型中的参数 若用多项式函数 可直接回归拟合 时间变量应从1开始 以免时间值过大 对模型产生不利影响 用指数函数时先取对数将其转化为线性函数再拟合 对搪瓷盆问题 若用皮尔曲线 可先估计出L 再转化为一个线性函数去估计参数 若L不易估计 可用三和值法 查找文献 去估计参数 当选择了多个模型 如何确定谁最适合 一看各个模型的拟合情况 即相关指标 如参数的显著性 拟合优度等 二是将拟合曲线与数据相比较 看它们的吻合程度 可作图直观观察 也可用平均相对误差 即式 要小于10才可接受 因为建模的目的主要是预测 近期数据的吻合情况更要关注 此外 也可预留部分数据拟合时不用 再用拟合所得模型比较这些数据与模型的预测值 此外 还要考虑序列的发展趋势以及模型的未来趋势 吻合者为好 有时 前期数据太多太久反而对拟合参数产生不利影响 也可考虑舍弃早期一些数据进行拟合 综合这些去选择模型 对例3 利用回归得到模型中的参数 从而将实际数据 曲折的线 与三个模型曲线画在同一个图中 上图 由图可以看出 三次曲线 虚线 与实际数据最接近 因而它最好 若计算MAPE 平均相对误差 比较后会得同一结论 但是 如果用极限等分析未来函数变化趋势 则三次曲线未必最佳 例5下表是我国1995 2000年各季度的商品零售总额 季节模型 Q 季度 对表中数据 做图如下 由图可知 销量既有线性增长趋势 又有明显的季节特征 故可以用模型yt a bt di i 1 2 3 4 其中di是季节增量 可由式di Di Di T Di m 1 T m得到 其中T是周期 m是年份数 Dt yt a bt 利用回归方法可得a bt 4901 46 164 88t 此时Dt yt a bt 的图像如下 由图知 Dt中已基本没有了增长趋势 再计算Dt的对应各季度的平均值 得 d1 d2 d3 d4 11 67 421 40 379 73 818 30 从而得趋势部分模型f t a bt di 对时间序列yt 当把趋势部分提炼出来 则其余下部分是yt f t xt 若它仅仅是由独立的随机干扰引起的误差 即xt N o 2 称为白噪声 则yt f t 就是我们要寻找的模型 可以用它进行预测 若xt 即模型误差不是白噪声 通常前后项间相关 则表明在其中仍有联系可循 需要把它们提炼出来 以便得到更好的模型 这就是下面的自回归模型 自回归 AR p 模型 时刻t的数值xt可表为前p个时刻的数值xt 1 xt 2 xt p的线性组合 再加上t时刻的白噪声 t 即xt a 1xt 1 2xt 2 pxt p t 9 其中a 通过序列减去其均值可化为0 因此可以没有 1 2 p是常数 模型说明xt的值与前面p期历史数据有关 注意与线性回归的差别 线性回归是不同变量间的关系 自回归是同一变量不同时期间的关系 另外 它还受随机扰动 t N 0 2 白噪声 各 t相互独立 并与前面的xt不相关 的影响 移动平均 MA q 模型 形式为xt t 1 t 1 2 t 2 q t q 10 6自回归模型 其中 1 2 q是常数 t t 1 t 2 t q是本期和前面若干期的随机扰动 模型含义 用一些时期的干扰的线性组合来表达当前值 自回归移动平均 ARMA p q 模型 xt 1xt 1 2xt 2 pxt p t 1 t 1 2 t 2 q t q 11 其中p和q分别是自回归和移动平均的阶数 模型含义 用前一些时期的数值及其干扰的线性组合来表达当前值 AR p 和MA q 是其特殊情形 7自回归模型识别及参数确定 用自回归建模要求序列 Xt 是 宽 平稳时间序列 即需满足 1 E Xt 与t无关 2 Var Xt 2 3 Cov Xt Xt k rk 与k有关 与t无关 由1 2 知 一个平稳时间序列 在图形上应表现出围绕均值的不断波动过程 因为随机 而非平稳时间序列则不然 往往表现出在不同时段有不同的均值 如持续上升或下降或有周期性 即有趋势或周期波动 7 1自相关和偏相关函数 自协方差函数rk E xt xt k 其中 是xt的均值 自相关函数 k rk r0 它反映的是xt与xt k的相关性 由概率知识知 k 1 该值的绝对值越大 说明两者的相关性越强 等于或接近于0说明没有相关性 实际计算中rk和 k的计算采用如下的表达式 n为样本容量 偏相关函数 是指在给定了Xt 1 Xt 2 Xt k 1的情况下Xt与Xt k的相关性 它反映的是在其他滞后期1 2 k 1的X已知的情况下Xt与Xt k的条件相关性 一般用符号 kk表示 并也有 kk 1 该值的绝对值越大 相关性越强 计算公式为 该式是由右端的矩阵方程得到 其中 k j k 1 j kk k 1 k j k 1 2 m j 1 2 k 1 此式是一个递推式子 虽然复杂 但在应用中k一般小于等于3 此外 有专门的命令语句可用 大家不必为计算而烦恼 7 2序列平稳性的检验 方法1图像观察平稳时间序列在图形上表现为围绕均值的不断波动过程 而非平稳时间序列则不然 方法2看自相关函数 k 或偏相关函数 kk 的图像平稳序列的 k是迅速衰减到0的 拖尾或截尾 非平稳序列的 k则是不衰减 缓慢衰减或 T T是季节周期 非常大 方法3根据得到的自回归模型中的式子xt 1xt 1 2xt 2 pxt p 视其为差分方程 若对应特征方程的根的绝对值都小于1 平稳 否则 不平稳 其他方法 略 注 k或 kk的拖尾和截尾 截尾 k在k q时全为0的性质称为q步截尾性 若它不能在某步之后截尾 而是随着k的增大而迅速衰减到0 受一负指数函数 如y e kx 控制 或如正弦函数似的震荡 称为拖尾性 此外 由于随机性 k全为0是不可能的 因此 截尾是指 k突然变的很小 并很接近于0 注 AR p 模型平稳的充要条件是它的p个特征根都在单位圆内 MA q 模型总是平稳的 ARMA p q 的平稳性与其AR p 部分相同 当序列 Xt 非平稳时 说明趋势f t 存在 除了季节趋势和明显的指数增长或阻滞增长趋势外 在短期内一般可用多项式函数近似 当为多项式函数时 如yt a bt t 通过不断的差分就可得到一个平稳序列 因此 对序列进行差分 季节规律用季节差分 是将非平稳序列变为一个平稳序列的常用手段 一般不超过两次 续例5 yt的图像明显上升 Dt具有明显周期性 都非平稳 它们的自相关系数图像如下 由图像可知 非平稳 图中有两条对称的蓝线 是随机变量的2 线 落在线内说明可以接受相关系数为0 95 的置信度 线外则不可 7 3白噪声检验 对于时间序列 yt 需要把其中的规律或项间关联全部提炼出来 使得残差 t 余下的部分 仅为一个白噪声 因此 检验残差序列是否为白噪声是判断模型是否合理以及建模是否需要终止的一个条件 设模型的残差序列为 t 记 计算 其中n为数据个数 m为最大时滞 m视数据多少取 n 4 n 10 或 n0 5 Qm近似服从 2 m 分布 对给定显著性水平 若Qm大于 2 m 则拒绝假设 否则接受 认为是白噪声 此外 t 是否为白噪声也可通过其相关函数来判断 若其 k和 kk都很小 可认为是 否则 不是 另 数模书中346页的投资问题给出了一个残差序列自相关性诊断方法 可画出 t t 1 r 1 n 1 r 2 n 的图像观察 也可用DW检验 仅检验一阶相关性 但一般够了 续例3 f t 用三次多项式拟合 22 372 f t 142 2 102 7t 7 67t2 0 22t3 残差序列为xt yt f t 其图像如右上 可以看出 序列基本平稳 但从 t t 1图像 右下 看 残差相邻项间有很强的相关性 用DW检验 算得的DW 0 9757 应通不过检验 此时Q8 22 372 也不满足白噪声检验 这说明残差序列间尚有信息留待提取 7 4自回归模型识别和参数确定 自回归有AR MA和ARMA三种模型 可从 k与 k的特性来判别 AR p 模型 k拖尾 kk滞后p阶后截尾 MA q 模型 k 滞后q阶后截尾 kk 拖尾 ARMA p q 模型 k 拖尾 kk 拖尾 续例5 将Dt进行季节差分 zt Dt 4 Dt 画图如下 由图像看 虽然zt没有明显的上升或下降趋势 以及季节特征 但并不很好的表现出围绕均值的波动 将zt再差分 仍记为zt 其图像如下 可以看出 序列平稳 检验自己做 了 下图是最后得到的zt的自相关函数和偏相关函数图像 由图知 自相关函数可以认为是拖尾的 偏相关函数则是截尾的 在2或4处 zt接近于白噪声 大家去检验 故应选AR模型 下图是例3的残差xt的自相关函数和偏相关函数图像 由图知 自相关函数可以认为是拖尾的 偏相关函数则是截尾的 在2处 故应选AR模型 显然不能认为是白噪声 k与 kk的截尾处的 严格 判断 k 若在某个q0 含 之前 k显著不为0 当q q0 q0 1 q0 2 q0 M中满足式的个数少于M的68 3 或上面不等式右端乘2 但比例变为95 5 则可近似认为 k在q0处截尾 其中N为数据个数 M同上 对 kk 判断方法类似 只是不等式是 kk 1 N0 5 或 kk 2 N0 5 截尾值q0可用来判别序列自回归或移动的阶数 若 k 要特别关注 与 kk既不拖尾也不截尾 说明序列非平稳 或有季节特征 需进行相关处理 虽然ARMA p q 模型具有一般性 但它也最复杂 另外 用ARMA p q 模型时 各 t k通常是未知的 不可观测量 因此 当用该模型时 必须求出前面的各 t 这不容易 对AR p 模型 其偏相关函数截尾的值 基本 就是回归的阶数 此外 也可用不同阶模式进行回归 残差平方和最小的值就是回归的阶数 看后面 对MA q 方法是类似的 对ARMA p q 模型 不能直接从相关函数得到大致的阶数 但残差平方和规则仍适用 方法是从低阶开始 向高阶拟合 在拟合的模型中选残差小者 或者遇到第一个残差可认为是白噪声的模式即停止 注意 当用高阶自回归移动平均模型去拟合序列时 拟合的效果总会提高的 不可能降低 即残差平方和会下降 但到了一定阶数后 阶数的再提高产生的效果会是微小的 非实质性的 这时的拟合属于过拟合 拟合过度 在建模中对模型还有一个 简约性 要求 即在精度相近的模型中我们要选择简单模型 为此 又有一个一般的定阶准则 AIC准则 记 a2 模型的剩余平方和 实际观察数据个数 模型中参数个数 则AIC p q log a2 2 p q n p q的确定应使AIC p q 达到最小 该式子应该也适合AR和MA模型 但这个方法必须对多种模型求参数 拟合 算残差 计算量大 故仅对ARMA用它 且尽量避免用它 7 3模型的定阶 确定了模型的阶数 就要确定其中的各个系数 一个常用的准则是残差平方最小准则 模型中各个系数的确定应使得用模型计算各个时刻的x值时 残差 实际值与计算值的差 的平方和达到最小 对AR p 模型 根据回归原理 可将xt 1 xt 2 xt p作为自变量 xt作为函数 用命令regress去做 还可用如下式子计算回归系数 对MA q 模型 一阶的仍可用回归 高阶的就不行了 只好根据定义来 比如对MA 2 有xt t 1 t 1 2 t 2 即 t xt 1 t 1 2 t 2 设给定序列 x1 x2 x15 我们令x1 1 则 2 x2 1 1 x2 1x1 3 x3 1 2 2 1 x3 1 x2 1x1 2x1 如此可得到各个 t 1 2的确定应使得下式最小 7 4回归系数的计算 此外 对MA q 模型 高阶的可用下式确定其系数这是非线性方程 求解也不容易 此外 专门软件如SAS中也许有命令计算 也许新版的matlab也有 大家查查看 对ARMA p q 模型 可用类似于MA q 的方法确定系数 当然 更复杂 也可用如下近似方法 先用回归方法求出自回归的系数 再用xt 1xt 1 2xt 2 pxt p t 1 t 1 2 t 2 q t q这个移动平均模型来做 当p q很大时模型会很复杂 计算也困难 但实际中这种模型是不多的 一般小于等于3 甚或2 除非能显著的减少误差 我们都尽量用简单模型来做 8模型预测及相关说明 得到一个模型 当然要用模型 而其中一个重要的应用就是预测 对AR模型xt 1xt 1