时间序列分析模型构建与MATLAB实现

式 （２）称 为 ｎ阶 自 回 归 模 型 ，记 为 ＡＲ（ｎ） 模型。 当 φi ＝０时 ，模 型 （１）变 为 ： ＝ at － θ1at −1 － θ 2 at − 2 － … － θ m at − m （ ３ ）
式 （３）称 为 ｍ阶 滑 动 平 均 模 型 ，记 为 ＭＡ （ｍ）模 型 。
的 是 ＡＩＣ定 阶 方 法 。 首 先 假 定 ＡＲＭＡ（ｐ，ｑ） 模 型 的 一 组 阶 数 ｋ，ｊ，然 后 利 用 ＡＲＭＡ（ｐ，ｑ） 模型的自回归逼近法求得白噪声的估计方 差 δ ２ ， 再 利 用 ＡＲＭＡ（ｐ，ｑ）模 型 的 ＡＩＣ定 阶 方 法 ，计 算 ＡＩＣ函 数 。 Ａ Ｉ Ｃ （ ｋ ， ｊ ） ＝ ｌ ｎ （δ ２ （ ｋ ， ｊ ） ） ＋ ２ （ ｋ ＋ ｊ ） ／ Ｎ ＡＩＣ（ｋ，ｊ）的 最 小 值 点 （ p ' ， q ' ）称 为 （ｐ，ｑ） 的 ＡＩＣ定 阶 。 也 就 是 ＡＲＭＡ（ｐ，ｑ）模 型 的 阶 数。 模型的参数估计： 求 模 型 参 数 的 方 法 很 多 ，有 矩 估 计 法 、 自回归逼近法、 最大似然估计法等。 这里我 们 讲 的 是 ＡＲＭＡ（ｐ，ｑ）模 型 的 自 回 归 方 法 。 对 零 均 值 化 后 的 观 测 数 据 x1 ， x2 ，… x N ｎ 拟 合 ＡＲＭＡ（ｐ，ｑ）模 型 时 ， 可 采 用 如 下 的 自 回归逼近方法。 首 先 为 数 据 建 立 ＡＲ模 型 。 曲子回归阶数的上届Ｐ＝［
X T X X T ε X TY β = T T T ε Y ε X ε ε
再估计出白噪声的方差：
−1
科技资讯 Ｓ Ｃ Ｉ Ｅ Ｎ Ｃ Ｅ ＆ Ｔ Ｅ Ｃ Ｈ Ｎ Ｏ Ｌ Ｏ Ｇ Ｙ Ｉ Ｎ Ｆ Ｏ Ｒ Ｍ Ａ Ｔ Ｉ Ｏ Ｎ
φi （ ｉ ＝ １ ， ２ ， … ， ｎ ） 自 回 归 参
xi
的 自 回 归 滑 动 平 均 模 型 ，记 为 ＡＲＭＡ（ｎ，ｍ）
xt ＝ φ1 xt −1 ＋ φ2 xt −2 ＋…＋ φn xt − n ＋ at （２）
参考文献
［１］ 何 书 元 ．应 用 时 间 序 列 分 析 ［Ｍ］．北 京 大 学 出 版 社 ，２００３． ［２］ 田 铮 ［译 ］．时 间 序 列 的 理 论 与 方 法 （第 ２ 版 ）［Ｍ］．高 等 教 育 出 版 社 ，施 普 林 格 出 版 社 ，２００１． 附录： ％假 设 ｐ１＝－０．１３，ｐ２＝２．０５，ｑ１＝１．３０， ｑ２＝－１．０８ ％时 间 序 列 模 型 为 ：Ｘｔ＝－０．１３（Ｘｔ－１） ＋２．０５（Ｘｔ－２）－１．３０（Ａｔ－１）＋１．０８（Ａｔ－２）＋Ａｔ Ｎｎ＝ｉｎｐｕｔ（＇请 输 入 观 测 值 个 数 Ｎｎ＇）； ｆｏｒ ｏ＝１：１０ ％本 循 环 用 于 统 计 ｐ１＝－０．１３； ｐ２＝２．０５； ｑ１＝１．３０； ｑ２＝－１．０８；
Q ( a , b) =
t = L +1
∑ (x − ∑ a x
t j =1
N
p
j t− j
− ∑ b j ε t' − j )
j =1
q
2
σ' =
2
Q ( a, b ) N −L
３ 模型的定阶、参数估计的 ＭＡＴＬＡＢ 实 现
程 序 见 附 录 ，程 序 开 始 通 过 产 生 服 从 正 态 分 布 的 随 机 序 列 ，并 假 设 模 型 参 数 ，以 便检验通过程序计算出的模型参数是否正 确。 运 行 程 序 之 后 ，得 出 的 结 果 为 ： ｎｎａ ＝ ２ ｎｎｂ ＝ ２ ｐｍａ ＝ －０．１３４６ ２．０５１５ ｐｍｂ ＝ １．３０２２ －１．０７３２ 证明程序正确。
θ1at −1 － θ 2 at − 2 － … － θ m at − m ＋ at （ １ ）
式 中 ，称 为
xt ＝ φ1 xt −1 ＋ φ2 xt −2 ＋ … ＋ φn xt − n －
４ 结语
通 过 对 时 间 序 列 模 型 用 ＭＡＴＬＡＢ的 编 程 ，明 白 了 时 间 序 列 模 型 的 作 用 及 应 用 范 围 ，更 加 请 除 了 时 间 序 列 的 建 模 过 程 ，时 间 序列的建模过程实际上是一个循环的过 程 ，也 是 一 个 统 计 的 过 程 ，对 以 后 相 关 数 据 的处理提供了一种有用的方法。
ａ＝ｒａｎｄ（１，Ｎｎ）； ｘ（１）＝ａ（１）； ｘ（２）＝ｐ１＊ｘ（１）－ｑ１＊ａ（１）＋ａ（２）； ｘ（３）＝ｐ１＊ｘ（２）＋ｐ２＊ｘ（１）－ｑ１＊ａ（２）－ｑ２＊ａ （１）＋ａ（３）； ｆｏｒ ｉ＝４：Ｎｎ ｘ（ｉ）＝ｐ１＊ｘ（ｉ－１）＋ｐ２＊ｘ（ｉ－２）＋ｐ３＊ｘ（ｉ－ ３）－ｑ１＊ａ（ｉ－１）－ｑ２＊ａ（ｉ－２）－ｑ３＊ａ（ｉ－３）＋ａ（ｉ）； ｅｎｄ ｘ ； ％＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ％利 用 信 息 量 原 则 定 阶 ／／／／／／ ％用 ＡＲ模 型 、 ｘ估 计 参 数 ｐ ｆｏｒ ｎ＝１：ｓｑｒｔ（Ｎｎ） ｘ２＝ｘ（（ｎ＋１）：Ｎｎ）＇； ｆｏｒ ｊ＝１：Ｎｎ－ｎ ｆｏｒ ｊｉ＝１：ｎ ｘ１（ｊ，ｊｉ）＝ｘ（ｎ－ｊｉ＋ｊ）； ｅｎｄ ｅｎｄ ｐ＝（ｉｎｖ（ｘ１（１：Ｎｎ－ｎ，：）＇＊ｘ１（１：Ｎｎ－ ｎ，：）））＊ｘ１（１：Ｎｎ－ｎ，：）＇＊ｘ２（１：Ｎｎ－ｎ）； ｐ； ｆｏｒ ｊ＝ｎ＋１：Ｎｎ ｆｏｒ ｉｉ＝１：ｎ ａｘ（ｉｉ）＝ｐ（ｉｉ）＊ｘ（ｊ－ｉｉ）； ｅｎｄ ｃｔ（ｊ）＝ｘ（ｊ）－ｓｕｍ（ａｘ）； ｅｎｄ ｂ（ｎ）＝ｓｕｍ（ｃｔ（ｎ＋１：Ｎｎ）．＾２）／（Ｎｎ－ｎ）； ＡＩＣ（ｎ）＝ｌｏｇ（ｂ（ｎ））＋２＊ｎ／Ｎｎ； ｅｎｄ ＡＩＣ； ％＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ％求 信 息 量 的 最 小 值 ｐｐ＝ＡＩＣ（１）； ｎｎ＝１； ｆｏｒ ｉ＝１：ｓｑｒｔ（Ｎｎ） ｉｆ ｐｐ＞ＡＩＣ（ｉ） ｐｐ＝ＡＩＣ（ｉ）； ｎｎ＝ｉ； ｅｎｄ ｅｎｄ ｎｎ ％＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ％求残差 ｆｏｒ ｊ＝ｎｎ＋１：Ｎｎ ｆｏｒ ｉｉ＝１：ｎｎ ａｘ（ｉｉ）＝ｐ（ｉｉ）＊ｘ（ｊ－ｉｉ）； ｅｎｄ ｃｔ（ｊ）＝ｘ（ｊ）－ｓｕｍ（ａｘ）； ｅｎｄ ａｔ０（ｎｎ＋１：Ｎｎ）＝ｃｔ（ｎｎ＋１：Ｎｎ）； ％估 计 参 数 ａａ，ｂｂ ｆｏｒ ｐ＝１：５ ｆｏｒ ｑ＝１：５ ｈｈ＝［ｐ，ｑ，ｎｎ］；
xL xL −1 Λ x x Λ X = L +1 L Μ Μ x N −1 x N − 2 Λ x L− p +1 xL− p + 2 Μ xN − p
数 ； θ j （ｊ＝１，２，… ，ｍ）称 为 滑 动 平 均 参 数 ； ｛ at ｝之 一 序 列 为 白 噪 声 序 列 ，（１）式 称 为 模型。 特 殊 地 ，当 θ j ＝０时 ，模 型 （１）变 为 ：
１ 时间序列基本思想概述
无论是按时间序列排列的观测数据还 是 按 空 间 位 置 顺 序 排 列 的 观 测 数 据 ，数 据 之间都或多或少的存在统计自相关现象。 时 间 序 列 分 析 是 ２０世 纪 ２０年 代 后 期 开 始 出 现 的 一 种 数 据 处 理 方 法 ，是 系 统 辨 识 与 系 统 分 析 的 重 要 方 法 之 一 ，是 一 种 动 态 的 数 据处理方法。 时 间 序 列 分 析 的 特 点 在 于：逐 次 的 观 测 值 通 常 是 不 独 立 的 ，且 分 析 必 须 考 虑 到 观 测 资 料 的 时 间 顺 序 ，当 逐 次 观 测 值 相 关 时 ，未 来 数 值 可 以 由 过 去 观 测 资 料 来 预 测 ，可 以 利 用 观 测 数 据 之 间 的 自 相 关 性建立相应的数学模型来描述客观现象的 动态特征。 时 间 序 列 分 析 的 基 本 思 想 是 ：对 于 平 稳、 正态、 零均值的时间序列 {xt } ，若 {xt } 的 取 值 不 仅 与 其 前 ｎ 步 的 各 个 取 值 {xt −1} ，
{xt −2 } ，… ， {xt − n } 有 关 ，而 且 还 与 前 ｍ步 的 干
扰 {at −1} ， {at − 2 } ， …
{at − m } 有 关 （ ｎ ， ｍ ＝ １ ，
２，… ），则 按 多 元 现 行 回 归 的 思 想 ， 可 得 到 最 一 般 的 ＡＲＭＡ模 型 ：
２ ＡＲＭＡ（ｎ，ｍ）模型的构建过程
ＡＲＭＡ（ｎ，ｍ） 模 型 建 立 的 一 般 步 骤 包 括 ：数 据 的 获 取 与 预 处 理 、 模型结构的选 择、 模型的定阶、 模型参数的估计、 模型检 验和模型预测等几个步骤。 下面重点叙述一下模型的定阶、 参数 估计的原理和方法。 模型的定阶： 模 型 定 阶 的 方 法 很 多 ，这 里 我 们 介ＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ
科技资讯
时间序列分析模型构建与 Ｍ Ａ Ｔ Ｌ Ａ Ｂ 实现
谭 （重庆市勘测院 重庆 ４０００２０） 摘 要 ：主要讲述了时间序列分析的基本思想，模型构建方法、 过程及步骤，并通过假设的时间序列用ＭＡＴＬＡＢ软件，实现了模型的定阶、 模型参数的估计过程。 关键词 ：时间序列分析 时间序列模型 ＭＡＴＬＡＢ方法实现 中 图 分 类 号 ：ＴＰ３ 文 献 标 识 码 ：Ａ 文 章 编 号 ：１６７２－３７９１（２００９）０９（ｂ）－０２５３－０２
２５３
科技资讯
２００９ ＮＯ．２６ ＳＣＩＥＮＣＥ ＆ ＴＥＣＨＮＯＬＯＧＹ ＩＮＦＯＲＭＡＴＩＯＮ
学 术 论 坛
Ｌ＝ｍａｘ（ｈｈ）； ｆｏｒ ｊ＝１：Ｎｎ－Ｌ ｆｏｒ ｋｍ＝１：ｐ ｘｘ１１（ｊ，ｋｍ）＝ｘ（Ｌ－ｋｍ＋ｊ）； ｐ； ｅｎｄ ｅｎｄ ｆｏｒ ｊ＝１：Ｎｎ－Ｌ ｆｏｒ ｊｍ＝１：ｑ ｅ１（ｊ，ｊｍ）＝ａｔ０（Ｌ－ｊｍ＋ｊ）； ｅｎｄ ｅｎｄ ａａｂｂ＝ｉｎｖ（［ｘｘ１１（１：Ｎｎ－Ｌ，：）＇＊ｘｘ１１（１： Ｎｎ－Ｌ，：） ｘｘ１１（１：Ｎｎ－Ｌ，：）＇＊ｅ１（１：Ｎｎ－ Ｌ，：）；ｅ１（１：Ｎｎ－Ｌ，：）＇＊ｘｘ１１（１：Ｎｎ－Ｌ，：） ｅ１ （１：Ｎｎ－Ｌ，：）＇＊ｅ１（１：Ｎｎ－Ｌ，：）］）＊［ｘｘ１１（１： Ｎｎ－Ｌ，：）＇＊ｘ（（Ｌ＋１）：Ｎｎ）＇；ｅ１（１：Ｎｎ－Ｌ，：） ＇＊ｘ（（Ｌ＋１）：Ｎｎ）＇］； ｐａｂ＝ａａｂｂ＇； ｐａ＝ｐａｂ（１：ｐ）； ｐｂ＝ｐａｂ（ｐ＋１：ｐ＋ｑ）； ｐａ； ｐｂ； ％＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ％计算信息量 ｆｏｒ ｊ＝Ｌ＋１：Ｎｎ ｆｏｒ ｉｉ＝１：ｐ ａｘ（ｉｉ）＝ｐａ（ｉｉ）＊ｘ（ｊ－ｉｉ）； ｅｎｄ ｆｏｒ ｋｋ＝１：ｑ ｂｘ（ｋｋ）＝ｐｂ（ｋｋ）＊ｘ（ｊ－ｋｋ）； ｅｎｄ ａｂｔｔ（ｊ）＝ｘ（ｊ）－ｓｕｍ（ａｘ）－ｓｕｍ（ｂｘ）； ｅｎｄ ａｂ（ｐ，ｑ）＝ｓｕｍ（ａｂｔｔ（Ｌ＋１：Ｎｎ）．＾２）／（Ｎｎ－ Ｌ）； ＡＩＣＢ（ｐ，ｑ）＝ｌｏｇ（ａｂ（ｐ，ｑ））＋２＊（ｐ＋ｑ）／Ｎｎ； ＡＩＣＢ； ｅｎｄ ｅｎｄ ％＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ％求 信 息 量 的 最 小 值 ｐｐｐｐ＝ＡＩＣＢ（１，１）； ｎｎａ＝１； ｎｎｂ＝１； ｆｏｒ ｉｋ＝１：５ ｆｏｒ ｊｋ＝１：５ ｉｆ ｐｐｐｐ＞ＡＩＣＢ（ｉｋ，ｊｋ） ｐｐｐｐ＝ＡＩＣＢ（ｉｋ，ｊｋ）； ｎｎａ＝ｉｋ； ｎｎｂ＝ｊｋ； ｅｎｄ ｅｎｄ ｅｎｄ ｎｎａ ｎｎｂ ％＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊ ％估 计 参 数 ａａ，ｂｂ ｐ＝ｎｎａ； ｑ＝ｎｎｂ； ｈｈ＝［ｐ，ｑ，ｎｎ］； Ｌ＝ｍａｘ（ｈｈ）； ｆｏｒ ｊ＝１：Ｎｎ－Ｌ ｆｏｒ ｋｍ＝１：ｐ ｘｘ１１（ｊ，ｋｍ）＝ｘ（Ｌ－ｋｍ＋ｊ）； ｐ； ｅｎｄ ｅｎｄ ｆｏｒ ｊ＝１：Ｎｎ－Ｌ ｆｏｒ ｊｍ＝１：ｑ ｅ１（ｊ，ｊｍ）＝ａｔ０（Ｌ－ｊｍ＋ｊ）； ｅｎｄ ｅｎｄ ａｂｃ＝ｉｎｖ（［ｘｘ１１（１：Ｎｎ－Ｌ，：）＇＊ｘｘ１１（１： Ｎｎ－Ｌ，：） ｘｘ１１（１：Ｎｎ－Ｌ，：）＇＊ｅ１（１：Ｎｎ－ Ｌ，：）；ｅ１（１：Ｎｎ－Ｌ，：）＇＊ｘｘ１１（１：Ｎｎ－Ｌ，：） ｅ１ （１：Ｎｎ－Ｌ，：）＇＊ｅ１（１：Ｎｎ－Ｌ，：）］）＊［ｘｘ１１（１： Ｎｎ－Ｌ，：）＇＊ｘ（（Ｌ＋１）：Ｎｎ）＇；ｅ１（１：Ｎｎ－Ｌ，：） ＇＊ｘ（（Ｌ＋１）：Ｎｎ）＇］； ｐｍｍ＝ａｂｃ＇； ｐｍａ＝ｐｍｍ（１：ｐ）； ｐｍｂ＝ｐｍｍ（ｐ＋１：ｐ＋ｑ）； ｐｍａ ｐｍｂ ｅｎｄ