时间序列 · 学习笔记一

明·煜

已于 2022-09-15 09:38:57 修改

阅读量1.7k

点赞数 3

分类专栏：数学建模学习笔记文章标签：学习算法机器学习

于 2022-09-07 21:08:20 首次发布

本文链接：https://blog.csdn.net/weixin_53195427/article/details/126753790

版权

数学建模学习笔记专栏收录该内容

7 篇文章 9 订阅

订阅专栏

时间序列

时间序列预测法
将预测对象按照时间顺序排列起来，构成一个所谓的时间序列，从所构成的这一组时间序列过来的变化规律，推断今后变化的可能性及变化趋势、变化规律，就是时间序列预测法
时间序列模型
时间序列模型其实是一种回归模型，其基于的基本原理是，一方面承认事物发展的延续性，运用过去时间序列的数据进行统计分析就能推测事物发展趋势；另一方面又充分考虑到偶然因素影响而产生的随机性,为了消除随机波动的影响，利用历史数据，进行统计分析，并对数据进行适当的处理，进行趋势预测。
对时间序列模型的评价

优点
能够充分运用原时间序列的各项数据，计算速度快，对模型参数有动态确定的能力，精度较好，采用组合的时间序列或者把时间序列和其他模型组合效果更好
缺点
不能反映事物的内在联系，不能分析两个因素的相关关系，只适用于短期预测

一、确定性时间序列分析方法

时间序列预测技术就是通过对预测目标自身时间序列的处理，来研究其变化趋势。一个时间序列往往是以下几类变化形式的叠加或耦合：

长期趋势变动：指时间序列朝着一定的方向持续上升或下降，或停留在某一水平上的倾向，他反映了客观事物的主要变化趋势。
季节变动
循环变动。通常是指周期为一年以上，由非季节因素引起的涨落起伏波形相近的波动
不规则变动。通常分为突变变动和随机变动

通常用 $T_t$ 表示长期趋势项， $S_t$ 表示季节变动趋势项， $C_t$ 表示循环变动趋势项， $R_t$ 表示随机干扰项。常用的确定性时间序列模型有以下几种类型

加法模型： $y_t = T_t+S_t+C_t+R_t$
乘法模型： $y_t = T_t \cdot S_t \cdot C_t \cdot R_t$
s
混合模型： $y_t = T_t \cdot S_t + R_t$ $y_t = S_t+T_t \cdot C_t \cdot R_t$ 式中： $y_t$ 为观测目标的观测记录，均值 $E(R_t)=0$ ,方差 $Var(R_t) = \sigma^2$

如果在预测时间范围以内，无突然变动且随即变动的方差 $\sigma^2$ 较小，并且有理由认为过去和现在的演变趋势将继续发展到未来时，可用一些经验方法进行预测。具体方法如下：

1. 移动平均法

1. 原理分析
二次移动平均值计算公式为 $M_t^{(2)} = \frac{1}{N}(M_t^{(1)}+\dots+M_{t-N+1}^{(1)})=M_{t-1}^{(2)}+\frac{1}{N}(M_t^{(1)}-M_{t-N}^{(1)})$
当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均方法建立预测模型： $\hat{y}_{t+1} = M_t^(1) = \frac{1}{N}(y_t+y_{t-1}+\dots+y_{t-N+1}),t = N,N+1,\dots,T$
其预测标准误差为： $\sqrt{\frac{\sum_{t=N+1}^T(\hat{y}-y_t)^2}{T-N}}$

那么，N值应该如何选择呢？

最近N期序列值的平均值作为未来各期的预测结果。一般N取值范围： $\leq N \leq 200$

当历史序列的基本趋势变化不大其序列中随机变动成分较多时，N值应该取大一些，否则N的取值应小一些。
在有确定的季节变动周期的资料中，移动平均的项数应取周期长度
选择最佳N值得一个有效方法时：比较若干模型的预测误差，预测标准误差最小值为好。

2. 案例分析

问题描述

某企业1~11月的销售收入时间序列如表所示，试用一次简单移动平均法预测12月的销售收入

月份 $t$	1	2	3	4	5	6	7	8	9	10	11	12
销售收入 $y_t$	533.8	574.6	606.9	649.8	705.1	772.0	816.4	892.7	963.9	1015.1	1102.7	?

分别取 $N = 4, N = 5$ 的预测公式
$\hat{y}_{t+1}^{(1)}=\frac{y_t+y_{t-1}+y_{t-2}+y_{t-3}}{4},t = 4,5,\dots,11$
$\hat{y}_{t+1}^{(2)}=\frac{y_t+y_{t-1}+y_{t-2}+y_{t-3}}{4},t = 4,5,\dots,11$
最终带入 $t = 11$ ,计算得到 $\hat{y}_{12}$ 的值以及相应的预测标准误差，发现 $N = 4$ 时的标准误差小于 $N = 5$ 时的标准误差，因此取 $N = 4$ 时的预测结果

MATLAB代码

  ```md
  y = [533.8 574.6 606.9 649.8 705.1 772.0 816.4 892.7 963.9 1015.1 1102.7];
  m = length(y);
  n = 4:1:11; % n为移动平均的项数
  n_best = 0;p = 100000; % 初始值选取一个无限大的值
  y_p = 0;
  for i = 1:length(n) % 循环计算两个N值情况下的预测值和预测标准误差,由于n值不同，用元胞数组yhat来记录数值
      for j = 1:m-n(i)+1
      yhat{i}(j) = sum(y(j:j+n(i)-1))/n(i); % 计算模型的各点的估计值
      end
      y12(i) = yhat{i}(end); % 提出12月的预测值
      s(i) = sqrt(mean((y(n(i)+1:end)-yhat{i}(1:end-1)).^2));% 求预测的标准误差
      if s(i)<p
      y_p = y12(i);
      n_best = n(i);
      p = s(i);
      end
  end
  y_p,n_best % 输出预测标准误差最小的预测值和所选N值的大小
  ```

3. 方法评价
移动平均法只适合做近期预测，而且是预测目标的发展取值变化不大的情况。如果目标的发展趋势存在其他变化，采用简单移动平均法就会产生较大的预测偏差和滞后

2. 指数平滑法

1. 综述

指数平滑法解决的问题：
一次移动平均实际上认为最近N期数据对未来的影响相同，都加权1/N,而N期以前的数据对未来值没有影响,加权为0。但是，二次及更高次移动平均的权数却不是1/N,且次数越高，权数的结构越复杂。一般来说，历史数据对未来值的影响是随时间间隔的增长而递减的。
因此，更切合实际的方法应该是对各期观测值依时间顺序进行加权平均作为预测值，指数平滑法可以满足这一条件，而且具有简单的递推形式。
指数平滑法的种类
指数平滑法根据平滑次数的不同，又分为一次指数平滑法、二次指数平滑法和三次指数平滑法

2. 一次指数平滑法

预测模型

在这里插入图片描述

指数平滑的本质在于加权系数符合指数规律，并且具有平滑数据的功能。
观察公式： $S_t^{(1)} = \alpha \sum_{j = 0}^{\infty}(1-\alpha)^jy_{t-j}$ 这表明 $S_t^{(1)}$ 是全部历史数据的加权和，加权系数分别为 $\alpha,\alpha(1-\alpha),\alpha(1-\alpha)^2,\dots$ ，因此加权系数符合指数规律
因此预测模型为： $\hat{y}_{t+1}=S_t^{(1)}$ $\hat{y}_{t+1 }=ay_t+(1-\alpha)\hat{y}_t$

加权系数的确定
在进行指数平滑时，加权系数的选择很重要，由递归公式可知：
$\hat{y}_{t+1 }=ay_t+(1-\alpha)\hat{y}_t$ $\alpha$ 的大小其实规定了在新预测值中新数据和原预测值所占的比重。可以把上式改写为新的式子 $\hat{y}_{t+1 } = \hat{y}_t+\alpha(y_t - \hat{y}_t)$ 因此可以得到，新预测值是根据预测误差对原预测值进行修正而得到的， $\alpha$ 值越大，修正幅度就越大； $\alpha$ 值越小，修正幅度就越小。

若选取 $\alpha = 0$ ,则 $\hat{y}_{t+1} = \hat{y}_t$ ,即下期预测值等于本期预测值，在预测过程中不考虑任何新信息；若选取 $\alpha = 1$ ,则 $\hat{y}_{t+1} = y_t$ ,即下期预测值就等于本期观测值，完全不相信过去的信息。

因此 $\alpha$ 应该根据时间序列的具体性质在0-1之间选择，并遵循以下原则：

如果时间序列波动不大，比较平稳（说明过去信息参考价值很大）， $\alpha$ 的值应该取小一点，如0.1~0.5。以减少修正幅度，使得预测模型能够包含较长时间序列的信息
如果时间序列具有迅速且明显的变动倾向，则 $\alpha$ 的取值应该大点，如0.6~0.8，使得预测模型灵敏度高一点，以便迅速跟上数据的变化。
对不同数据的选取，我们选择预测标准误差最小。其计算公式为： $\sqrt{E[(y-\hat{y})^2]}$ 即根号下误差均值的平方
1. 初始值确定
  用一次指数平滑法进行预测，需要选择合适的初始值 $s_0^{(1)}$ 。初始值是由预测者估计或指定的。但时间序列的数据较多狮，比如在20个以上时，初始值对以后的预测值的影响很小，可以选择第一期数据为初始值。如果时间序列的数据较小时（在20个以内），初始值对以后的预测值影响很大，这时就必须认证研究何时正确确定初始值。一般以最初几期实际值的平均值作为初始值
4.案例分析：
1976年-1987年某种电器的销售额如表所示，试预测1988年该电器的销售额：

解题思路如下：

显然在取不同的 $\alpha$ 值是，得到预测值是不尽相同的。因此需要比较他们的预测标准误差，最终得到的结果是：

代码如下：

yt = xlsread('data1.xlsx'); % 实际销售额数据以列向量的形式储存在.xlsx文件中
n = length(yt);alpha = [0.2 0.5 0.8];m = length(alpha) % 获得权重向量
yhat(1,[1:m]) = (yt(1)+yt(2))/2; % 根据公式确定初始值
for i = 2:n
    yhat(i,:) = alpha*yt(i-1)+(1-alpha).*yhat(i-1,:);
end
err = sqrt(mean((repmat(yt,1,m)-yhat).^2)); % 计算预测标准误差
xlswrite('output1',yhat);
yhat1988 = alpha*yt(n)+(1-alpha).*yhat(n,:)

3. 二次指数平滑法

但当时间序列的变动出现直线趋势时，用一次指数平滑法进行预测，仍存在明显的滞后偏差。因此必须加以修正，建立二次指数平滑，利用之后偏差的规律建立直线趋势模型
在这里插入图片描述

3. 差分指数平滑法

略

4. 具有季节性特点的时间序列的预测

使用季节系数法，具体流程如下：
案例分析

代码如下：

data2 = xlsread('data2.xlsx'); % 读入数据，行标为年份，列表为季度
a = data2; % 与分析所用的字母一致
[m,n] = size(a);
a_mean = mean(mean(a)); % 计算所有数据的算术平均值
aj_mean = mean(a); % 计算同季度的平均值
bj = aj_mean/a_mean; % 计算季度系数
w = 1:m; % w_i为第i年的权数，按自然数列取值
yhat = w*sum(a,2)/sum(2); % 预测下一年的年加权平均值。
% 此处的sum（a，2）是对a进行行求和，得到的是列向量。跟w这个行向量相乘后得到一个1*1矩阵
yjmean = yhat/n; % 计算预测年份的季平均值
yjhat = yjmean*bj; % 预测年份的季节预测值

二、平稳时间序列模型

1. 时间序列的基本概念

$A RM A$ 时间序列模型是一种重要的平稳时间序列模型。 $A RM A$ 时间序列分为三种类型：
$A R$ 序列，即自回归序列
$M A$ 序列，即移动平均序列
$A RM A$ 序列，即自回归移动平均序列

AR ( p ) 序列

AR模型总结：
引入后移算子 $B$ 时，记算子多项式 $\varphi(B) = 1 - \varphi_1B-\varphi_2B^2-\dotsb-\varphi_pB^p$ 该式也可写为： $\varphi(B)X_t = \varepsilon_t$

MA(q)序列

MA模型总结
一个平稳白噪声通过一个线性系统，如果其输出是白噪声的线性叠加那么这一输出服从MA模型
引入后移算子 $B$ 时，记算子多项式 $\varphi(B) = 1 - \theta_1B-\theta_2B^2-\dotsb-\theta_pB^p$ 该式也可写为： $X_t = \theta(B)\varepsilon_t$
ARMA(p,q)序列

2. ARMA模型的构建及预报

在实际的建模问题中，分为以下步骤：

进行模型的识别与定阶，即要判断AR§,MA(q)或ARMA(p,q)模型的类别，并要估计阶数p，q。
当模型定阶完成后，就要对模型参数 $\pmb{\varphi} = [\varphi_1,\varphi_2,\dots,\varphi_n]^T和\pmb{\theta} = [\theta_1,\theta_2,\dots,\theta_n]^T$ 进行估计
模型参数估计完后，需要对模型进行检验，即要检验 $\varepsilon$ 是否是平稳白噪声。
模型建立完成后，最主要的应用是进行预报

（1）. ARMA模型的构建

ARMA模型定阶AIC准则
ARMA模型的参数估计
MATLAB的内置函数
ARMA模型检验的 $\chi^2$ 检验

（2）. ARMA(p,q)序列的预报

在这里插入图片描述

三、时间序列的MATLAB相关工具箱的命令

1. 系统辨识工具箱相关命令的使用

案例一：求解AR(2)模型,并利用模型进行预测

在这里插入图片描述

代码如下

clc,clear
elps = randn(10000,1);x(1:2) = 0; % 初始化第一个数据和第二个数据为0
for i = 3:10000
    x(i) = -0.6*x(i-1)-0.2*x(i-2)+elps(i); % 产生模拟数据
end
x = x'; m =ar(x,2) % 进行参数估计，之前以及利用Matlab软件识别得到AR(2)模型
xhat = forecast(m,x,3) % 计算3个预测值

MATLAB的计算结果显示为：
在这里插入图片描述

因此得到的预测结果是：
在这里插入图片描述

案例二：通过枚举法利用AIC准则来对模型定则
通过枚举法利用AIC准则，来对模型定阶，选取AIC取值最小的p，q
作为模型的阶数。

在这里插入图片描述


clc;clear;
randn('state',sum(clock)); % 初始化随机数发生器
elps = randn(1,10000); % 产生10000个服从标准正态分布的随机数
x(1) = 0; % 赋初始值
for j = 2:10000
    x(j) = 0.8*x(j-1)+elps(j)-0.4*elps(j-1); % 产生样本点
end
x = x'; % 转换成列向量
for i = 0:3
    for  j = 0:3
        if i==0 & j==0
            continue % arma(p,q)模型中，p，q不能同为0
        end
        m = armax(x,[i,j]); % 拟合模型，已知数据必须是列向量
        myaic = aic(m); % 计算AIC指标
        fprintf('p = %d,q = %d,AIC = %f\n',i,j,myaic); % 显示计算结果
    end
end

% 以上程序可以帮助我们穷举所有AIC值，找到最小的AIC值时的p，q输入模型
p = input('输入阶数p = '); q = input('输出阶数q = '); % 输入模型的阶数
m = armax(x,[p,q]);  % 拟合指定参数的p，q的模型
xp = predict(m,x); % 计算已知数据的预测值
res = x - xp; % 计算残差向量
h = lbqtest(res) % 进行chi2检验（即Ljung-Box检验）
    
% lbqtest检验结果是，h = 0时通过检验

四、ARMA序列与季节性序列

前面讲述了平稳时间序列的建模方法，在实际中遇到的时间序列却不一定是平稳的，往往有3个特性，即趋势性，季节性，平稳性。需要采用差分方法（即Box - Jenkins方法）消除其趋势性和季节性，最终得到平稳序列，并假设为ARMA模型。

1. ARIMA 序列及其预报

通过差分方法可以使序列变得平稳

在这里插入图片描述

ARIMA(p,q,d)序列

在这里插入图片描述

那么应该如何确定自己的模型已经平稳了呢？
在这里插入图片描述

那么问题又来了，拖尾和截尾又是什么意思呢
目前我并没有找到数理的证明方法，来说明某个序列是拖尾的或者截尾的，只能够通过观察图像来进行分析。即绘制ACF和PACF图
具体的图像可以看这篇文章：https://blog.csdn.net/Caiqiudan/article/details/118059325

ARIMA的预报

在这里插入图片描述

2. 案例分析

在这里插入图片描述

补充资料 AIC准则和BIC准则

在这里插入图片描述

AIC和BIC都是越小越好，共同点是都利用了最大似然估计

（1）计算过程

在这里插入图片描述

关于差分次数的确定，即为什么d = 1
可以观察下r21和r22的图像，两个图像模样差不多，下面只画出一种：

在这里插入图片描述

证明了自相关函数和偏相关函数是截尾的

（2）MATLAB代码

clc,clear;
%% 模型建立
a = textread('hua.txt'); % 把原始数据按照原来的排列格式存放在纯文本文件hua.txt。a的值并不是排列成一个列向量，而是一个矩阵
a = nonzeros(a'); % 按照原来的顺序去掉0元素;v = nonzeros(A)返回A中非零元素的满列向量，v中元素按列排列

% 这里取a'的原因在于，由于v =nonzeros(A)，会按照列的顺序输出非0元素。
% 即先判断输出的应是[A(1,1),A(2,1),A(3,1),...,A(M,N)]^T
% 而根据题意知，表格的时间先后是一行一行来的，因此先要对a进行装置

r11 = autocorr(a); % 计算自相关函数
r12 = parcorr(a); % 计算偏相关函数
da = diff(a); % 计算一阶差分
r21 = autocorr(da); % 计算自相关函数
r22 = autocorr(da) % 计算偏相关函数
% 通过计算相关系数和偏相关系数，确定ARIMA（p,q,d）中d = 1；
% 下面绘制自相关函数和偏相关函数的图像，从而观察其是否截尾或拖尾
figure(1);stem(r21);
figure(2);stem(r22);

n = length(da); % 计算差分后的数据个数
k = 0; % 初始化试探模型的个数
for i = 0:3
    for j = 0:3
        if i ==0 & j == 0
            continue % ARIMA(p,q,d)中的p和q不能同时为0
        elseif i ==0
            ToEstMd = arima('MALags',1:j,'Constant',0); % 指定模型的结构，当i=0时，变为MA模型 
        elseif j ==0
            ToEstMd = arima('ARLags',1:i,'Constant',0); % 指定模型结构，当j=0时，变成AR模型
        else
            ToEstMd = arima('ARLags',1:i,'MALags',1:j,'Constant',0); % 指定模型结构
        end
        k = k+1; % 试验模型数+1
        R(k) = i;M(k) = j;
        [EstMd,EstParamCov,logL,info] = estimate(ToEstMd,da); % 对模型进行拟合，即对参数进行估计
        numParams = sum(any(EstParamCov)); % any 测试每一列是否有非零元素,如果有则返回1；
        % 根据Akaike和贝叶斯信息进行定阶（即使AIC和BIC）
        [aic(k),bic(k)] = aicbic(logL,numParams,n); % 确定第k个模型的AIC和BIC信息
    end
end
fprintf('R,M,AIC,BIC的对应值如下\n % f'); % 显示计算结果
check = [R',M',aic',bic'];

% 通过上述模型可以辨别穷举所有的情况，我们需要选择AIC和BIC最小的来作为我们的模型。

%% 模型的应用-预测
r = input('输入阶数R = '); m = input('输出阶数M = ');
ToEstMd = arima('ARLags',1:r,'MALags',1:m,'Constant',0); % 指定模型结构
dx_Forecast = forecast(EstMd,10,'Y0',da); % 计算10步的预测值
x_Forecast = a(end)+cumsum(dx_Forecast); %

% 如果A是一个向量， cumsum(A) 返回一个向量，该向量中第m行的元素是A中第1行到第m行的所有元素累加和；
% 这是因为dx_Forecast获得得是前一项相当于后一样的增加值，因此我们才需要在x_Forecast中做这样的处理

五、季节性序列及其预报

由季节性因素或其他周期性因素引起的周期性变化的时间序列，称为季节性时间序列，相应的模型为季节性模型。
一般地，对周期 $s$ 的序列，可先进行差分运算，即 $\nabla_sX_t = (1-B^s)X_t$ $\nabla_s^dX_t = (1-B^s)^dX_t$

1. 案例分析

在这里插入图片描述

解答
在这里插入图片描述

MATLAB代码

clc;clear
%% 模型建立
x = load('water.txt'); % 原始数据将会按照表中的格式存放在纯文本文件water.txt中
x= x';x=x(:); % 按照时间得到先后次序，把数据变成列向量
s = 12; % 周期 s =12
n =12; % 预报数据个数
m1 = length(x); % 原始数据的个数
for i = s+1:m1
    y(i-s) = x(i)-x(i-s); % 进行周期差分变换
end
w = diff(y); % 消除趋势性的差分运算。
% 注意，此时y的元素个数已经与x的元素个数不同了，应该比其少12；w的元素个数也与y不同，应该比y少1
% 下面绘制自相关函数和偏相关函数的图像，从而观察其是否截尾或拖尾
r21 = autocorr(w); % 计算自相关函数
r22 = autocorr(w) % 计算偏相关函数
figure(1);stem(r21);
figure(2);stem(r22);
m2 = length(w); % 计算最终差分后数据的个数
k = 0; % 初始化试探模型的个数
for i = 0:3
    for j = 0:3
        if i ==0 & j==0
            continue % 因为ARIMA(p,q,d)中的pq不能同时为0
        elseif i ==0 
            ToEstMd = arima('MAlags',1:j,'Constant',0); % 指定模型的结构
        elseif j==0
            ToEstMd = arima('ARlags',1:i,'Constant',0); % 指定模型的结构
        end
        k = k+1;R(k) = i;M(k) = j;
        [EstMd,EstParamCov,logL,info] = estimate(ToEstMd,w'); % 对模型进行拟合，即对参数进行估计
        numParams = sum(any(EstParamCov)); % any 测试每一列是否有非零元素,如果有则返回1；
        % 根据Akaike和贝叶斯信息进行定阶（即使AIC和BIC）
        [aic(k),bic(k)] = aicbic(logL,numParams,n); % 确定第k个模型的AIC和BIC信息
    end
end
fprintf('R,M,AIC,BIC的对应值如下\n % f'); % 显示计算结果
check = [R',M',aic',bic'];

% 通过上述模型可以辨别穷举所有的情况，我们需要选择AIC和BIC最小的来作为我们的模型。

%% 模型预测
r = input('输入阶数R = '); m = input('输出阶数M = ');
ToEstMd = arima('ARLags',1:r,'MALags',1:m,'Constant',0); % 指定模型结构
w_Forecast = forecast(EstMd,n,'Y0',w'); % 计算12步的预测值
yhat = y(end)+cumsum(w_Forecast); % 求一阶差分的还原值

% 由于还做了周期性差分，因此还需要进行一次还原
for j = 1:n
    x(m1+j) = yhat(j)+x(m1+j-s);
end
xhat = x(m1+1:end) % 截取n个预测值

关于w的自相关函数和偏相关函数如下图
在这里插入图片描述

明·煜

关注

3
点赞
踩
17

收藏

觉得还不错? 一键收藏
2
评论
时间序列 · 学习笔记一

时间序列预测法将预测对象按照时间顺序排列起来，构成一个所谓的时间序列，从所构成的这一组时间序列过来的变化规律，推断今后变化的可能性及变化趋势、变化规律，就是时间序列模型，其基于的基本原理是，一方面承认事物发展的延续性，运用过去时间序列的数据进行统计分析就能推测事物发展趋势；另一方面又充分考虑到偶然因素影响而产生的随机性,为了消除随机波动的影响，利用历史数据，进行统计分析，并对数据进行适当的处理，进行趋势预测。对时间序列模型的评价优点。
复制链接

扫一扫