时间序列预测之ARMAARIMA序列及季节性序列matlab实现

Posted 2022-02-12 Icy Hunter

ARMA是一种平稳时间序列模型，即均值和协方差不随时间的平移而改变。
ARMA有三种类型

AR序列

MA序列

ARMA序列

但是由于ARMA只能处理平稳序列，而现实中的问题往往有趋势性或周期性等。为了得到平稳序列，我们对数据进行差分运算，使得新序列成为平稳序列，就能够进行ARMA分析，因此ARIMA模型，是在ARMA的基础上多了差分运算，使得其能够处理的序列范围增加了。
ARIMA序列


例题1：

clc,clear
a = [17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 17.0
    16.7 17.4 17.2 17.4 17.4 17.0 17.3 17.2 17.4 16.8
    17.1 17.4 17.4 17.5 17.4 17.6 17.4 17.3 17.0 17.8
    17.5 18.1 17.5 17.4 17.4 17.1 17.6 17.7 17.4 17.8
    17.6 17.5 16.5 17.8 17.3 17.3 17.1 17.4 16.9 17.3
    17.6 16.9 16.7 16.8 16.8 17.2 16.8 17.6 17.2 16.6
    17.1 16.9 16.6 18.0 17.2 17.3 17.0 16.9 17.3 16.8
    17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 0
    17.3 17.4 17.7 16.8 16.9 17.0 16.9 17.0 16.6 16.7
    16.8 16.7 16.4 16.5 16.4 16.6 16.5 16.7 16.4 16.4
    16.2 16.4 16.3 16.4 17.0 16.9 17.1 17.1 16.7 16.9
    16.5 17.2 16.4 17.0 17.0 16.7 16.2 16.6 16.9 16.5
    16.6 16.6 17.0 17.1 17.1 16.7 16.8 16.3 16.6 16.8
    16.9 17.1 16.8 17.0 17.2 17.3 17.2 17.3 17.2 17.2
    17.5 16.9 16.9 16.9 17.0 16.5 16.7 16.8 16.7 16.7
    16.6 16.5 17.0 16.7 16.7 16.9 17.4 17.1 17.0 16.8
    17.2 17.2 17.4 17.2 16.9 16.8 17.0 17.4 17.2 17.2
    17.1 17.1 17.1 17.4 17.2 16.9 16.9 17.0 16.7 16.9
    17.3 17.8 17.8 17.6 17.5 17.0 16.9 17.1 17.2 17.4
    17.5 17.9 17.0 17.0 17.0 17.2 17.3 17.4 17.4 17.0
    18.0 18.2 17.6 17.8 17.7 17.2 17.4 0 0 0]';

a = nonzeros(a);
r11 = autocorr(a);
r12 = parcorr(a);
da = diff(a);
r21 = autocorr(da)
r22 = parcorr(da)
n = length(da);
k = 0;
for i = 0:3
    for j = 0:3
        if i == 0 & j == 0
            continue
        elseif i == 0
            ToEstMd = arima('MALags', 1:j, 'Constant', 0);
        elseif j == 0
            ToEstMd = arima('ARLags', 1:i, 'Constant', 0);
        else 
            ToEstMd = arima('ARLags', 1:i, 'MALags', 1:j, 'Constant', 0);
        end
        k = k + 1;
        R(k) = i;
        M(k) = j;
        [EstMd, EstParamCov, logL, info] = estimate(ToEstMd, da);
        numParams = sum(any(EstParamCov));
        [aic(k), bic(k)] = aicbic(logL, numParams, n);
    end
end
fprintf('R, M, AIC, BIC的对应值如下\\n ');
check = [R', M', aic', bic']
r = input('输入阶数R=');
m = input('输入阶数M=');
ToEstMd = arima('ARLags', 1:r, 'MALags', 1:m, 'Constant', 0);
[EstMd, EstParamCov, logL, info] = estimate(ToEstMd, da);
dx_Forecast = forecast(EstMd, 10, 'Y0', da)
x_Frecast = a(end) + cumsum(dx_Forecast)

我代码敲出来结果有一点偏差，可能数据打错了几个吧，太多了，不想看了。

例题2

clc,clear
x = [9.40 8.81 8.65 10.01 11.07 11.54 12.73 12.43 11.64 11.39 11.1 10.85
    10.71 10.24 8.48 9.88 10.31 10.53 9.55 6.51 7.75 7.8 5.96 5.21
    6.39 6.38 6.51 7.14 7.26 8.49 9.39 9.71 9.65 9.26 8.84 8.29
    7.21 6.93 7.21 7.82 8.57 9.59 8.77 8.61 8.94 8.4 8.35 7.95
    7.66 7.68 7.85 8.53 9.38 10.09 10.59 10.83 10.49 9.21 8.66 8.39
    8.27 8.14 8.71 10.43 11.47 11.73 11.61 11.93 11.55 11.35 11.11 10.49
    10.16 9.96 10.47 11.70 10.1 10.37 12.47 11.91 10.83 10.64 10.29 10.34];
x = x';
x = x(:);
s = 12; % 周期
n = 12; % 预测数
m1 = length(x); % 原始数据数
for i = s+1:m1
    y(i-s) = x(i) - x(i-s); % 进行周期差分变换
end
w = diff(y); % 消除趋势性差分运算
m2 = length(w); % 计算最终差分后数据个数
k = 0; % 模型个数

for i = 0:3
    for j = 0:3
        if i == 0 & j == 0
            continue
        elseif i == 0
            ToEstMd = arima('MALags', 1:j, 'Constant', 0);
        elseif j == 0
            ToEstMd = arima('ARLags', 1:i, 'Constant', 0);
        else 
            ToEstMd = arima('ARLags', 1:i, 'MALags', 1:j, 'Constant', 0);
        end
        k = k + 1;
        R(k) = i;
        M(k) = j;
        [EstMd, EstParamCov, logL, info] = estimate(ToEstMd, w'); % 模型拟合
        numParams = sum(any(EstParamCov)); % 计算拟合参数个数
        [aic(k), bic(k)] = aicbic(logL, numParams, m2);
    end
end
fprintf('R, M, AIC, BIC的对应值如下\\n ');
check = [R', M', aic', bic']
r = input('输入阶数R=');
m = input('输入阶数M=');
ToEstMd = arima('ARLags', 1:r, 'MALags', 1:m, 'Constant', 0); % 指定模型结构
[EstMd, EstParamCov, logL, info] = estimate(ToEstMd, w'); % 拟合模型
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);% 求x的预测值
end
xhat = x(m1 + 1: end) % 截取n个预报值

分析：

这题我和他的结果是对的上的。

使用ARIMA模型简单来说就是先差分使得序列成为平稳序列，然后就是遍历AR和MA的参数的事了，看aic和bic越小的，就选哪个参数就好了。

本文参考的是司守奎，孙兆亮主编的数学建模算法与应用（第二版）