机械学习基础-10.从时间序列数据中学习-数据建模与机械智能课程自留

本节目标：

• 理解如何使用普通最小二乘法（OLS）将时间序列模型拟合到数据中。
• 知晓如何使用预测未来时间序列数据。
• 能够对自回归模型（AR）、移动平均模型（MA）、自回归 - 外生变量模型（ARX）和自回归移动平均 - 外生变量模型（ARMAX）进行分类。
• 知晓如何迭代地拟合自回归移动平均 - 外生变量模型（ARMAX）。

时间序列数据（TIME SERIES DATA）

1. 静态特征数据：之前在监督学习问题中考虑的是“静态”特征数据。输入是 $\hat{x}\in {\mathbb{R}}^m$，通过模型 $f(\hat{x})$ 得到预测值 $\hat{y}$。这表明输入是一个 $m$ 维的向量，经过模型计算后产生预测结果。
2. 动态（时间序列）特征数据：接下来考虑的是“动态”（按时间排序）的特征数据。输入是一个由 $N$ 个时间点的数据组成的向量 $\left[\begin{array}{c}\hat{x}(t)\\ ⋮\\ \hat{x}(t+N-1)\end{array}\right]\in {\mathbb{R}}^N$，经过模型 $f\left(\left[\begin{array}{c}\hat{x}(t)\\ ⋮\\ \hat{x}(t+N-1)\end{array}\right]\right)$ 得到预测值 $\hat{y}(N)$。这里强调了输入数据是按时间顺序排列的一系列数据点，并且特征数据包含“已知”量，如历史标签数据或外生输入。

动态系统术语

基本元素

• 输入/外生变量（Input/Exogenous variable）：表示为 $u(t)$，是系统的外部输入。
• 系统（System）：表示为 $x(t)$，是整个动态系统的核心部分。
• 干扰/误差（Disturbance/error）：表示为 $ϵ(t)$，是影响系统的外部干扰因素或误差项。
• 输出（Output）：表示为 $y(t)$，是系统产生的可测量结果。系统的内部状态通常用 $x(t)$ 表示，输出通常是内部状态的函数，即 $y(t)=g(x(t))$。当 $y(t)=x(t)$ 时，这种情况被称为全状态反馈（full state feedback）。

系统模型类型

• 连续时间模型（Continuous time models）：形式为 $\dot{x}(t)=f(x(t),u(t))$，其中线性系统 $\dot{x}(t)=Ax(t)+Bu(t)$ 是其一种特殊情况。这里 $\dot{x}(t)$ 表示 $x(t)$ 对时间的导数。
• 离散时间模型（Discrete time models）：形式为 $x(t+1)=f(x(t),u(t))$，描述了系统在离散时间点上的状态转移关系。

自回归模型（AR）

线性预测模型

线性预测模型假设
假设未来输出服从线性参数模型：$$y(t)=\psi (t{)}^{\top }\theta +ϵ(t)$$其中：

• $\psi (t)\in {\mathbb{R}}^m$ 是在时间 $t$ 已知的数据向量。
• $\theta \in {\mathbb{R}}^m$ 是我们可以进行优化的模型参数。
• $ϵ(t)$ 是建模和测量误差。

线性模型示例
$$y(t)={\theta }_{1}u(t)+{\theta }_2$$
在这个例子中，对应的数据向量 $\psi (t)=\left(\begin{array}{c}u(t)\\ 1\end{array}\right)$ 。模型参数 $\theta$ 是线性的。

具体来说，$\psi (t{)}^{\top }\theta$ 是向量 $\psi (t)$ 与向量 $\theta$ 的内积形式，展开后是 $\theta$ 各个分量的线性组合 。

例如，若 $\psi (t)=[{\psi }_1(t),{\psi }_2(t),\cdots ,{\psi }_m(t){]}^{\top }$，$\theta =[{\theta }_1,{\theta }_2,\cdots ,{\theta }_m{]}^{\top }$，则 $\psi (t{)}^{\top }\theta ={\sum }_{i=1}^m{\psi }_i(t){\theta }_i$ ，没有出现 $\theta$ 的非线性运算（如平方、开方、指数等），所以该模型是线性模型。

预测

模型拟合
假设我们已使用数据拟合了如下模型：$\hat{y}(t)=\psi (t{)}^{\top }\hat{\theta }$ 。这里 $\hat{y}(t)$ 是预测输出，$\psi (t)$ 是已知的数据向量，$\hat{\theta }$ 是估计的模型参数。

步骤
我们可以用这个模型对未来若干时间步的输出进行预测。预测结果 $\hat{Y}$ 是一个向量，其形式为 $$\hat{Y}=\left[\begin{array}{c}\hat{y}(1)\\ ⋮\\ \hat{y}(N)\end{array}\right]=\left[\begin{array}{c}\psi (1{)}^{\top }\hat{\theta }\\ ⋮\\ \psi (N{)}^{\top }\hat{\theta }\end{array}\right]=\left[\begin{array}{ccc}{\psi }_1(1)& \cdots & {\psi }_m(1)\\ ⋮& \ddots & ⋮\\ {\psi }_1(N)& \cdots & {\psi }_m(N)\end{array}\right]\hat{\theta }$$ 在预测的情境下，右侧的矩阵被称为回归矩阵（regressor matrix）。

示例
还是上面那个例子：$y(t)=\widehat{{\theta }_1}u(t)+\widehat{{\theta }_2}$ ，对应的数据向量 $\psi (t)=\left(\begin{array}{c}u(t)\\ 1\end{array}\right)$ 。假设 $u(t)$ 在时间步 $t=1,2,3,4,5$ 的取值分别为 $0,1,2,3,4$ ，则预测矩阵 $\hat{Y}=\left[\begin{array}{cc}0& 1\\ 1& 1\\ 2& 1\\ 3& 1\\ 4& 1\end{array}\right]\hat{\theta }$ 。该示例展示了从给定模型到确定数据向量，再到构建预测矩阵的完整过程。

使用延迟数据进行预测

背景
在时间阶段 $t$，我们可能会记得之前的输出 $y(1),\dots ,y(t-1)$。因此，我们可能希望将这些数据作为预测模型的一部分。这种模型被称为自回归（Auto - Regressive, AR）模型。

公式
AR模型记为 $AR(n_y)$,其表达式为 $$y(t)={\theta }_{1}y(t-1)+\cdots +{\theta }_{n_y}y(t-n_y)+ϵ(t)$$

其中，$y(t)$ 是当前时刻的输出，$y(t-1),\dots ,y(t-n_y)$ 是过去时刻的输出，${\theta }_1,\dots ,{\theta }_{n_y}$ 是模型参数，$ϵ(t)$ 是误差项。

$AR(n_y)$ 模型可以写成 $\hat{y}(t)=\psi (t{)}^{\top }\hat{\theta }$ 的形式，其中 $\psi (t)=[y(t-1),\dots ,y(t-n_y){]}^{\top }\in {\mathbb{R}}^{n_y}$。这类模型被称为自回归模型，因为“标签”数据同时也是“特征”数据。

在实际应用中，通常存在数据延迟收集的情况，所以AR模型非常有用。

自回归模型的回归矩阵

1. AR模型表达式
   考虑一个 $AR(n_y)$ 模型，其表达式同上，为 $$\hat{y}(t)=\widehat{{\theta }_1}y(t-1)+\cdots +\widehat{{\theta }_{n_y}}y(t-n_y)$$

2. 数据假设
   通常假定对于所有 $t<0$ 的情况，数据向量 $\psi (t)=0$ ，这其中包括 $y(t)$ 和 $u(t)$ 。时间 - 输出值如下：

t	0	…	N - 1
y(t)	y(0)	…	y(N - 1)

3. 回归矩阵构建
   预测向量 $\hat{Y}=\left[\begin{array}{c}\hat{y}(1)\\ ⋮\\ \hat{y}(N)\end{array}\right]$ ，它可以表示为 $\hat{Y}=\Psi \hat{\theta }$ 的形式，其中 $\Psi$ 是回归矩阵。
   

回归矩阵 $\Psi$ 的构建方式如下：

• 矩阵的第一行是 $[y(0),0,\cdots ,0]$ ，表示预测 $\hat{y}(1)$ 时，只用到了 $y(0)$ 这一个历史值（假设 $n_y\ge 1$ ），其他滞后项为 $0$ （基于 $t<0$ 时数据为 $0$ 的假设）。
• 第二行是 $[y(1),y(0),\cdots ,0]$ ，表示预测 $\hat{y}(2)$ 时，用到了 $y(1)$ 和 $y(0)$ 两个历史值。
• 以此类推，最后一行是 $[y(N-1),y(N-2),\cdots ,y(N-n_y)]$ ，表示预测 $\hat{y}(N)$ 时，用到了从 $y(N-1)$ 到 $y(N-n_y)$ 共 $n_y$ 个历史值（如果 $N-n_y\ge 0$ ）。

自回归 - 外生变量模型（ARX）

ARX模型定义

我们可以将AR模型扩展为ARX模型，记为 $ARX(n_y,n_u)$，其表达式为 $$y(t)={\theta }_{1}y(t-1)+\cdots +{\theta }_{n_y}y(t-n_y)+w_{1}u(t-1)+\cdots +w_{n_u}u(t-n_u)+ϵ(t)$$

其中，$y(t)$ 是当前时刻的输出，$y(t-1),\cdots ,y(t-n_y)$ 是过去时刻的输出，${\theta }_1,\cdots ,{\theta }_{n_y}$ 是对应输出的模型参数；$u(t-1),\cdots ,u(t-n_u)$ 是过去时刻的外生输入变量，$w_1,\cdots ,w_{n_u}$ 是对应外生输入的模型参数；$ϵ(t)$ 是误差项。

超参数说明
$n_y$ 和 $n_u$ 是自然数，属于“手动选取”的超参数，而非通过优化得到。$n_y$ 和 $n_u$ 越大，意味着决策变量越多，这会使模型对训练数据的拟合效果更好，但同时也面临过拟合的风险。

模型形式转换
$ARX(n_y,n_u)$ 模型可以写成 $\hat{y}(t)=\psi (t{)}^{\top }\hat{\theta }$ 的形式，其中：

• $\psi (t)=[y(t-1),\cdots ,y(t-n_y),u(t-1),\cdots ,u(t-n_u){]}^{\top }\in {\mathbb{R}}^{n_y+n_u}$ ，是由过去输出值和过去外生输入值组成的数据向量。
• $\hat{\theta }=[\widehat{{\theta }_1},\cdots ,\widehat{{\theta }_{n_y}},\widehat{w_1},\cdots ,\widehat{w_{n_u}}{]}^{\top }\in {\mathbb{R}}^{n_y+n_u}$ ，是由对应参数组成的参数向量。

ARX模型的回归矩阵

• 数据假设：通常假定对于所有$t<0$的情况，数据向量$\psi (t)=0$，这其中包括$y(t)$和$u(t)$。从$t=0$到$t=N-1$时刻的$y(t)$和$u(t)$的取值情况如下：

t	0	…	N - 1
y(t)	y(0)	…	y(N - 1)
u(t)	u(0)	…	u(N - 1)

• 预测向量与回归矩阵关系：预测向量$\hat{Y}=\left[\begin{array}{c}\hat{y}(1)\\ ⋮\\ \hat{y}(N)\end{array}\right]$，它可以表示为$\hat{Y}=\Psi \hat{\theta }$的形式，其中$\Psi$是回归矩阵。需注意，$\hat{Y}$向量使用的时间数据是 { 1 , ⋯   , N } \{1,\cdots,N\} {1,⋯,N}，而回归矩阵$\Psi$使用的时间数据是 { 0 , ⋯   , N − 1 } \{0,\cdots,N - 1\} {0,⋯,N−1}。
• 回归矩阵构建：回归矩阵$\Psi$是一个$N\times (n_y+n_u)$的矩阵。
  
  其构建方式为：矩阵中前半部分是基于过去输出$y(t)$的数据，从第一行开始，每行依次是不同时刻往前追溯$n_y$个时刻的$y(t)$值（不足部分补0）；后半部分是基于过去外生输入$u(t)$的数据，每行依次是不同时刻往前追溯$n_u$个时刻的$u(t)$值（不足部分补0）。

使用普通最小二乘法（OLS）拟合ARX模型

拟合目标
假设我们想用给定数据拟合一个$ARX(n_y,n_u)$模型。为实现这一目标，需要最小化误差平方和，其数学表达式为$$\min \sum _{i=1}^N(y_i-{\hat{y}}_i{)}^2=\min \Vert Y-\hat{Y}{\Vert }_2^2=\min \Vert Y-\Psi \theta {\Vert }_2^2$$ 其中，$y_i$是实际输出值，${\hat{y}}_i$是模型预测输出值，$Y$是实际输出值向量，$\hat{Y}$是预测输出值向量，$\Psi$是回归矩阵，$\theta$是模型参数向量。

该问题本质上就是OLS问题，存在解析解。而最小二乘问题的解析解为$${\theta }^{\ast }=({\Psi }^{\top }\Psi {)}^{-1}{\Psi }^{\top }Y$$其中$\Psi$是上面给出的ARX模型的回归矩阵。

在使用最小二乘法拟合模型时外生输入（激励输入）的作用

一、 最小二乘法条件
为了使用最小二乘法的解析公式，要求${\Psi }^{\top }\Psi$ 是可逆的。这是使用最小二乘法求解模型参数的一个重要前提条件。

在ARX模型这类情境下，不用外生输入，${\Psi }^{\top }\Psi$不一定可逆，原因如下：

• 从ARX模型结构来看
  ARX（自回归 - 外生变量）模型包含了系统自身过去输出以及外生输入变量对当前输出的影响。当只有系统自身过去输出（类似AR模型情况，没有外生输入）时，回归矩阵$\Psi$只由过去输出值构成。如果系统本身的动态特性导致这些过去输出值之间存在较强的相关性（比如输出序列呈现出高度的自相关性，在一些平稳的时间序列中较为常见），那么${\Psi }^{\top }\Psi$的行列式可能为零或者非常接近于零，从而导致${\Psi }^{\top }\Psi$不可逆。

• 从数据生成角度
  外生输入的一个重要作用是激励系统产生多样化、不相关的数据。如果没有外生输入，系统的输出仅依赖于自身过去状态，可能会陷入某种周期性或相关性较强的模式，使得$\Psi$矩阵的列向量之间存在线性相关关系。而${\Psi }^{\top }\Psi$可逆的一个重要条件是$\Psi$的列向量线性无关。当外生输入存在时，通过合理设计（如采用阶跃函数、正弦函数和等作为外生输入），可以打破这种相关性，增加$\Psi$矩阵列向量的线性独立性，进而保证${\Psi }^{\top }\Psi$可逆，以便能够使用最小二乘法的解析公式来求解模型参数。

二、丰富模型信息
外生输入可以引入额外信息，除了考虑系统自身过去的输出（如 AR 模型），还能纳入其他外部因素对系统输出的影响，使模型更贴近实际情况，提升对系统行为的解释和预测能力。例如在预测销售额时，除了依据过去的销售数据，还可将促销活动、季节等外生因素作为外生输入纳入模型。

三、系统激励与数据生成
作为激励系统的手段，外生输入能够促使系统产生不同的响应，生成丰富多样的数据。这些数据可用于更好地训练和验证模型，提高模型的泛化能力和准确性。不同类型的外生输入（如阶跃函数、正弦函数和）能从不同角度激发系统特性，帮助挖掘系统的动态行为。

常见外生输入（激励输入）

• 阶跃函数：$u(t)=\{\begin{array}{ll}0& \text{if }t<0\\ u_0& \text{otherwise}\end{array}$，即在$t<0$ 时取值为$0$，其他时刻取值为$u_0$ 。
• 正弦函数和：$u(t)={\sum }_{j=1}^N{a}_j\sin ({\omega }_{j}t+{\nu }_j)$ ，由不同幅度$a_j$ 和权重的正弦函数叠加而成。

激励输入的选择通常被视为一个超参数，一般通过启发式方法来选择。

拟合ARX模型的数值示例

1. 设定真实模型
   假设真实模型满足以下$ARX(1,1)$模型：$y(t)=0.5y(t-1)+0.75u(t-1)+ϵ(t)$，其中$ϵ(t)\in N(0,1)$，即误差项$ϵ(t)$服从均值为$0$、方差为$1$的正态分布。

2. 生成数据集
   使用随机激励序列生成数据集，对应的Matlab代码为u = 2*(1 - 2*rand(N, 1))。这行代码通过rand(N, 1)生成一个$N$行$1$列的随机数矩阵（随机数取值在$0$到$1$之间），再经过变换生成取值在$-2$到$2$之间的随机激励序列$u$。

3. 生成输出数据
   考虑正态分布噪声，生成输出数据$y$。Matlab代码如下：

e = randn(N, 1);
y = zeros(N, 1);
for i = 2:N
    y(i) = 0.5*y(i - 1)+0.75*u(i - 1)+e(i);
end

上述代码首先通过randn(N, 1)生成一个$N$行$1$列的服从标准正态分布的噪声向量$e$；然后初始化输出向量$y$为全$0$的$N$行$1$列矩阵；最后通过循环，根据设定的$ARX(1,1)$模型生成输出数据$y$。

4. 构建回归矩阵
   构建回归矩阵$\Psi$，Matlab代码为for i = 2:N; psi(i, :) = [y(i - 1) u(i - 1)]; end。该循环从第$2$个时间步开始，将前一个时间步的输出$y(i-1)$和外生输入$u(i-1)$组合成一行，构建回归矩阵$\Psi$。

5. 求解OLS问题
   使用普通最小二乘法（OLS）求解模型参数$\theta$，Matlab代码为theta = inv(psi'*psi)*psi'*y。这行代码根据最小二乘法的解析解公式${\theta }^{\ast }=({\Psi }^{\top }\Psi {)}^{-1}{\Psi }^{\top }Y$，计算得到模型参数$\theta$，其中inv表示求逆矩阵，psi'表示对psi矩阵进行转置。

在实际中，随着时间增加，数据点数量增加，模型参数误差会减小。

横轴为数据点数量，纵轴为参数误差。从图中可以看出，随着数据点数量的增加，参数误差总体呈下降趋势。在数据点数量达到一定程度（不是很多数据点之后），误差接近零。这表明，更多的数据点有助于提高模型参数估计的准确性，使得估计参数更接近真实参数。

真实输出与估计输出的比较

图表的横轴为时间步（Time steps），纵轴为输出（Output）。图中有两条曲线，蓝色曲线代表估计输出（Estimated output），橙色曲线代表真实输出（True Output）。可以看到，两条曲线的走势大致相似，但存在一定的波动差异。

计算公式：

• 真实输出$y_{true}(t)=0.5y_{true}(t-1)+0.75u(t-1)+ϵ(t)$ ，其中包含了系统自身过去的输出、外生输入以及误差项$ϵ(t)$ 。
• 估计输出$\hat{y}(t)=0.5001y_{true}(t-1)+0.7506u(t-1)$ ，这里的参数$0.5001$和$0.7506$是通过普通最小二乘法（OLS）并使用200,000个数据点估计得到的。

由于测量误差的存在，真实输出和预测输出之间总是存在微小差异。

具有内部动态的示例

假设真实模型满足以下式子：

• $x(t)=0.5x(t-1)+0.75u(t-1)$ ，描述了内部状态变量 $x(t)$ 随时间的变化，它与前一时刻的内部状态 $x(t-1)$ 和外生输入 $u(t-1)$ 相关。
• $y(t)=x(t)+ϵ(t)$ ，其中 $ϵ(t)\sim N(0,1)$ ，表示输出 $y(t)$ 由内部状态 $x(t)$ 和服从均值为 $0$ 、方差为 $1$ 的正态分布的噪声 $ϵ(t)$ 组成。模型存在内部动态，并且能够将内部状态作为带有小噪声的输出进行测量。

该模型与之前的系统非常相似，通过重新整理可得：$$y(t)=0.5y(t-1)-0.5ϵ(t-1)+0.75u(t-1)+ϵ(t)$$
进一步尝试将其写成线性框架形式：$$y(t)=\psi (t{)}^{\top }\theta -0.5ϵ(t-1)+ϵ(t)$$

该模型与之前的数值示例不同，因为输出测量受到延迟噪声的影响。之前的示例可能没有这种延迟噪声对输出测量的干扰情况。

病态数值示例：普通最小二乘法失效的情况

尝试使用OLS拟合以下模型：$$y(t)=0.5y(t-1)+0.75u(t-1)+ϵ(t)$$

其中，$y(t)$ 是当前时刻的输出，$y(t-1)$ 是前一时刻的输出，$u(t-1)$ 是前一时刻的外生输入，$ϵ(t)$ 和 $ϵ(t-1)$ 是噪声项。预测模型形式为 $\hat{y}(t)=\psi (t{)}^{\top }\hat{\theta }$ ，其中 $\psi (t)=[y(t-1),u(t-1){]}^{\top }$ 。

在这个模型中忽略了延迟噪声 $-0.5ϵ(t-1)$ ，这可能是导致OLS失效的原因

生成数据的代码

• 首先通过 u=2*(1 - 2*rand(N, 1)) 生成随机的外生输入序列 $u$ 。
• 然后用 e=randn(N, 1) 生成服从标准正态分布的噪声序列 $e$ ，并初始化输出序列 $y$ 为全零向量 y=zeros(N, 1) 。
• 通过循环 for i = 2:N ，按照模型公式 y(i)=0.5*y(i - 1)+0.75*u(i - 1)+e(i)-0.5*e(i - 1) 生成输出数据 $y$ 。
• 再通过另一个循环 for i = 2:N 构建回归矩阵 psi ，最后用 theta=inv(psi'*psi)*psi'*y 尝试使用OLS求解模型参数 theta 。

输出图像

图表展示了使用200,000个数据点时模型的表现，横轴为时间步（Time steps），纵轴为输出（Output）。图中有两条曲线，蓝色曲线代表估计输出（Estimated output），橙色曲线代表真实输出（True Output）。可以看到，估计输出和真实输出之间存在较大偏差，说明在这种情况下OLS方法未能很好地拟合模型。

原因：参数估计不收敛到零
下图是随着估计参数$\theta$与输入数据的关系图。

上图的横轴为数据点数量，纵轴为参数误差。从图中可以看出，无论使用多少数据点，参数误差始终没有减小到0.2以下，即误差不会随着数据点数量的增加而无限减小趋近于零。

在使用200,000个数据点后，估计参数为$\hat{\theta }=[0.2501,0.7495{]}^{\top }$，而真实参数为${\theta }^{\ast }=[0.5,0.75{]}^{\top }$。可以发现，估计参数与真实参数并不接近，存在一定的偏差。

似乎如果我们忽略延迟噪声$0.5ϵ(t-1)$，并应用普通最小二乘法（OLS），那么在参数估计中总是会存在一些误差。这表明，延迟噪声的存在对OLS方法的参数估计产生了显著影响，导致估计结果无法很好地收敛到真实参数值，即使使用大量数据也不能消除这种误差。

移动平均模型（MA）

普通最小二乘法解析解中的偏差

真实模型设定
假设真实的潜在动态关系由$y(t)=\psi (t{)}^{\top }{\theta }^{\ast }+ϵ(t)$给出，其中：

• $\psi (t)\in {\mathbb{R}}^m$是在时刻$t$已知的数据向量。
• ${\theta }^{\ast }\in {\mathbb{R}}^m$是真实的模型参数。
• $ϵ(t)$是“白噪声”，具有零均值且与过去或未来时刻不相关。

该式可以写成矩阵形式：$Y=\Psi {\theta }^{\ast }+ϵ$ 。

OLS的无偏性引理
OLS无偏性的引理推导过程：

• $\mathbb{E}(\hat{\theta })=\mathbb{E}(({\Psi }^{\top }\Psi {)}^{-1}{\Psi }^{\top }Y)$，其中$\hat{\theta }$是OLS估计的参数。
• 将$Y=\Psi {\theta }^{\ast }+ϵ$代入上式，得到$\mathbb{E}(({\Psi }^{\top }\Psi {)}^{-1}{\Psi }^{\top }(\Psi {\theta }^{\ast }+ϵ))$。
• 经过展开和化简，因为$\mathbb{E}(ϵ)=0$，所以$\mathbb{E}(\hat{\theta })={\theta }^{\ast }+({\Psi }^{\top }\Psi {)}^{-1}{\Psi }^{\top }\mathbb{E}(ϵ)={\theta }^{\ast }$，这表明在一定条件下OLS估计是无偏的，即估计参数的期望值等于真实参数。

存在的问题
如果$\Psi$不包含任何噪声（可能仅依赖于外生输入），那么最小二乘法会给出无偏估计。但在具有延迟噪声的自回归模型中，这个问题尤其麻烦。因为在这种情况下，$\Psi$的元素会依赖于$ϵ(t-1)$，导致$\mathbb{E}(ϵ(t-1{)}^2)\ne 0$的项出现在OLS公式中，从而使OLS估计产生偏差，不再具有无偏性。

自回归 - 滑动平均 - 外生变量模型（ARMAX MODELS）

我们可以将AR模型推广到自回归 - 滑动平均 - 外生变量（ARMAX）模型。

ARMAX模型的形式
一个$ARMAX(n_y,n_{ϵ},n_u)$模型可以写成$y(t)=\psi (t{)}^{\top }\theta +ϵ(t)$的形式，其中：

• $\psi (t)=[y(t-1),\cdots ,y(t-n_y),ϵ(t-1),\cdots ,ϵ(t-n_{ϵ}),u(t-1),\cdots ,u(t-n_u){]}^{\top }\in {\mathbb{R}}^{n_y+n_{ϵ}+n_u}$，是由过去输出值、过去噪声值和过去外生输入值组成的数据向量。
• $\theta =[{\theta }_1,\cdots ,{\theta }_{n_y},{\alpha }_1,\cdots ,{\alpha }_{n_{ϵ}},w_1,\cdots ,w_{n_u}{]}^{\top }\in {\mathbb{R}}^{n_y+n_{ϵ}+n_u}$，是对应的模型参数向量。

滑动平均（MA）项对应于模型依赖于延迟误差/噪声这一事实。之所以称为MA，是出于历史原因。

对ARMAX模型应用OLS来拟合模型参数会导致有偏估计。因此，接下来需要研究一种替代方法来将ARMAX模型拟合到数据中。

用于拟合ARMAX模型的迭代普通最小二乘法

基本思想：通过递归估计数据中的噪声来迭代更新模型参数。

初始化
以有偏估计进行初始化：$${\hat{\theta }}_1:=({\Psi }_1^{\top }{\Psi }_1{)}^{-1}{\Psi }_1^{\top }Y$$
其中${\Psi }_1:=[{\hat{\psi }}_1(0),\cdots ,{\hat{\psi }}_1(N){]}^{\top }$，
${\hat{\psi }}_1(t):=[y(t-1),\cdots ,y(t-n_y),u(t-1),\cdots ,u(t-n_u){]}^{\top }$

迭代步骤
从$n=1$到$N_{stop}$进行以下操作：

1. 估计误差/噪声：使用${\hat{\theta }}_n$估计每个时间阶段的误差/噪声，公式为$${\hat{e}}_n(t):=y(t)-{\hat{y}}_n(t)=y(t)-{\hat{\psi }}_n(t{)}^{\top }{\hat{\theta }}_n$$ 即当前时刻的真实输出$y(t)$减去基于当前参数估计${\hat{\theta }}_n$和数据向量${\hat{\psi }}_n(t)$得到的预测输出${\hat{y}}_n(t)$ 。
2. 求解OLS：使用${\Psi }_{n+1}:=[{\hat{\psi }}_{n+1}(0),\cdots ,{\hat{\psi }}_{n+1}(N){]}^{\top }$求解OLS，其中$${\hat{\psi }}_{n+1}(t):=[y(t-1),\cdots ,y(t-n_y),{\hat{e}}_n(t-1),\cdots ,{\hat{e}}_n(t-n_e),u(t-1),\cdots ,u(t-n_u){]}^{\top }$$ 此时的数据向量${\hat{\psi }}_{n+1}(t)$相比初始化时加入了上一步估计得到的延迟噪声项${\hat{e}}_n(t-1),\cdots ,{\hat{e}}_n(t-n_e)$ 。
3. 更新参数估计：更新参数估计${\hat{\theta }}_{n+1}:=({\Psi }_{n+1}^{\top }{\Psi }_{n+1}{)}^{-1}{\Psi }_{n+1}^{\top }Y$ 。通过上述迭代过程，不断改进对ARMAX模型参数的估计。

使用OLS拟合ARMAX模型的MATLAB伪代码

设定ARMAX模型
要拟合的为$y(t)=0.5y(t-1)+0.75u(t-1)+ϵ(t)-0.5ϵ(t-1)$ 。

生成输入 - 输出数据

u=2*(1 - 2*rand(N, 1));
e=randn(N, 1);
y=zeros(N, 1);
for i = 2:N
    y(i)=0.5*y(i - 1)+0.75*u(i - 1)+e(i)-0.5*e(i - 1);
end

上述代码首先通过u=2*(1 - 2*rand(N, 1))生成取值在 - 2到2之间的随机外生输入序列$u$ ；然后用e=randn(N, 1)生成服从标准正态分布的噪声序列$e$ ，并初始化输出序列$y$为全零向量y=zeros(N, 1) ；最后通过循环，依据给定的ARMAX模型生成输出数据$y$ 。

初始噪声估计

e_est=y - psi*theta

根据模型参数的初始近似值来估计噪声。这里的psi和theta需要提前有初始设定，该步骤用于得到初始的噪声估计值e_est 。

递归估计系统参数和噪声

for j = 1:N_stop
    for i = 2:N
        psi(i, :)=[y(i - 1) u(i - 1) e_est(i - 1)];
    end
    theta=inv(psi'*psi)*psi'*y;
    e_est=y - psi*theta;
end

通过两层循环进行迭代：

• 内层循环构建数据矩阵psi ，它包含了前一时刻的输出y(i - 1) 、前一时刻的外生输入u(i - 1)以及前一时刻估计的噪声e_est(i - 1) 。
• 接着使用OLS公式theta=inv(psi'*psi)*psi'*y更新模型参数theta 。
• 最后根据更新后的参数重新估计噪声e_est=y - psi*theta 。通过这样的递归过程，不断改进对系统参数和噪声的估计。

拟合ARMAX模型的数值示例

上图的横轴为数据点数量，纵轴为参数误差。从图中可以看出，随着数据点数量的增加，参数误差总体呈下降趋势。

通过迭代求解OLS并估计噪声$ϵ(t)$ ，当迭代停止条件$N_{stop}=10$时，估计的模型参数收敛到真实模型参数${\theta }^{\ast }=[0.5,0.75,-0.5{]}^{\top }$ 。这表明通过这种迭代OLS的方法能够有效地对ARMAX模型进行拟合，随着数据量的积累和迭代的进行，参数估计越来越准确，最终趋近于真实参数值。

手写示例

问题设定
考虑一个ARMAX(1, 1, 0)模型，初始化参数${\theta }_0=[21{]}^{\top }$，初始噪声估计${\widehat{ϵ}}_0=[10.5-1{]}^{\top }$。任务是估计系统噪声/误差，并利用它进行一次普通最小二乘法（OLS）迭代来更新参数。

模型形式为$$\hat{y}(t)=ay(t-1)+b\widehat{ϵ}(t-1)$$其中$\hat{y}(t)$是预测输出，$y(t-1)$是前一时刻的实际输出，$\widehat{ϵ}(t-1)$是前一时刻的噪声估计，$a$和$b$是待估计参数。

数据如下：

t	0	1	2	3
y(t)	-1	2	3	2

系统噪声估计
利用${\theta }_0=[ab{]}^{\top }$估计系统噪声，噪声估计向量${\psi }_0(t{)}^{\top }=[y(t-1),{\widehat{ϵ}}_0(t-1){]}^{\top }$。计算过程如下：

• $\widehat{{ϵ}_1}(0)=y(0)-{\psi }_0(0{)}^{\top }{\theta }_0=-1-0=-1$
• $\widehat{{ϵ}_1}(1)=y(1)-{\psi }_0(1{)}^{\top }{\theta }_0=2-2\times (-1)-1=3$
• $\widehat{{ϵ}_1}(2)=y(2)-{\psi }_0(2{)}^{\top }{\theta }_0=3-2\times 2-0.5=-1.5$

回归矩阵计算
计算回归矩阵${\Psi }_1$：

$${\Psi }_1=\left[\begin{array}{cc}y(0)& {\widehat{ϵ}}_1(0)\\ y(1)& {\widehat{ϵ}}_1(1)\\ y(2)& {\widehat{ϵ}}_1(2)\end{array}\right]$$
$$=\left[\begin{array}{cc}-1& -1\\ 2& 3\\ 3& -1.5\end{array}\right]$$

计算参数更新

计算${\Psi }_1^{\top }{\Psi }_1$
$${\Psi }_1^{\top }{\Psi }_1=\left[\begin{array}{ccc}-1& 2& 3\\ -1& 3& -1.5\end{array}\right]\left[\begin{array}{cc}-1& -1\\ 2& 3\\ 3& -1.5\end{array}\right]=\left[\begin{array}{cc}14& 2.5\\ 2.5& 12.25\end{array}\right]$$
这里是将回归矩阵${\Psi }_1$与其转置${\Psi }_1^{\top }$相乘，得到一个$2\times 2$的矩阵。

计算$({\Psi }_1^{\top }{\Psi }_1{)}^{-1}$
$$({\Psi }_1^{\top }{\Psi }_1{)}^{-1}=\left[\begin{array}{cc}cc-1& -1\\ 0.0741& -0.0151\\ -0.0151& 0.0847\end{array}\right]$$
即对上一步得到的矩阵求逆，得到逆矩阵。蓝色框内文字提示“Note, we do not use y(0) in Y”，即注意在$Y$中不使用$y(0)$。

计算${\Psi }_1^{\top }Y$
$${\Psi }_1^{\top }Y=\left[\begin{array}{ccc}-1& 2& 3\\ -1& 3& -1.5\end{array}\right]\left[\begin{array}{c}2\\ 3\\ -2\end{array}\right]=\left[\begin{array}{c}-2\\ 10\end{array}\right]$$

这里$Y$是由输出值组成的向量（根据提示未包含$y(0)$），将${\Psi }_1^{\top }$与$Y$相乘得到一个$2\times 1$的向量。

计算更新后的参数${\theta }_1$
$$\begin{array}{rl}{\theta }_1& =({\Psi }_1^{\top }{\Psi }_1{)}^{-1}{\Psi }_1^{\top }Y=\left[\begin{array}{cc}0.0741& -0.0151\\ -0.0151& 0.0847\end{array}\right]\left[\begin{array}{c}-2\\ 10\end{array}\right]=\left[\begin{array}{c}-0.2995\\ 0.8775\end{array}\right]\end{array}$$
根据OLS的参数更新公式，将前面计算得到的$({\Psi }_1^{\top }{\Psi }_1{)}^{-1}$与${\Psi }_1^{\top }Y$相乘，得到更新后的参数向量${\theta }_1$。

不同预测模型总结

自回归模型（AR）
表达式为$$AR(n_y):y(t)={\theta }_{1}y(t-1)+\cdots +{\theta }_{n_y}y(t-n_y)+ϵ(t)$$ 该模型基于当前时刻的输出$y(t)$是其过去若干时刻$n_y$输出值的线性组合，再加上一个白噪声项$ϵ(t)$。

移动平均模型（MA）
表达式为$$MA(n_{ϵ}):y(t)={\theta }_1ϵ(t-1)+\cdots +{\theta }_{n_{ϵ}}ϵ(t-n_{ϵ})+ϵ(t)$$此模型中当前时刻的输出$y(t)$是过去若干时刻$n_{ϵ}$噪声值的线性组合，再加上当前时刻的噪声。

自回归 - 移动平均模型（ARMA）
表达式为$$ARMA(n_y,n_{ϵ}):y(t)={\theta }_{1}y(t-1)+\cdots +{\theta }_{n_y}y(t-n_y)+w_1ϵ(t-1)+\cdots +w_{n_{ϵ}}ϵ(t-n_{ϵ})+ϵ(t)$$它综合了AR和MA模型的特点，既考虑了过去输出值的影响，也考虑了过去噪声值的影响。

带外生变量的自回归模型（ARX）
表达式为$$ARX(n_y,n_u):y(t)={\theta }_{1}y(t-1)+\cdots +{\theta }_{n_y}y(t-n_y)+w_{1}u(t-1)+\cdots +w_{n_u}u(t-n_u)+ϵ(t)$$该模型在AR模型基础上，加入了外生变量$u(t)$，即考虑了外部输入对输出的影响。

带外生变量的自回归 - 移动平均模型（ARMAX）
表达式为$ARMAX(n_y,n_{ϵ},n_u):y(t)={\theta }_{1}y(t-1)+\cdots +{\theta }_{n_y}y(t-n_y)+{\alpha }_1ϵ(t-1)+\cdots +{\alpha }_{n_{ϵ}}ϵ(t-n_{ϵ})+w_{1}u(t-1)+\cdots +w_{n_u}u(t-n_u)+ϵ(t)$ 。它结合了ARMA和ARX模型的特性，既包含过去输出、过去噪声的影响，也考虑了外生变量的作用。

非线性自回归 - 移动平均模型（NARMAX）
表达式为$$NARMAX(n_y,n_{ϵ},n_u):y(t)=F(y(t-1),\cdots ,y(t-n_y),ϵ(t-1),\cdots ,ϵ(t-n_{ϵ}),u(t-1),\cdots ,u(t-n_u))$$与前面的线性模型不同，NARMAX模型是一个非线性模型，$F$表示一个非线性函数，它可以捕捉更复杂的动态关系。

超参数选择
应使用交叉验证（Cross - validation）来选择每个模型的超参数，以优化模型性能。超参数是在模型训练之前需要设定的参数，交叉验证通过将数据集划分为多个子集进行多次训练和验证，来评估不同超参数组合下模型的表现，从而选择最优的超参数。