第十讲-时间序列预测

习题8.4

1. 题目要求

2.解题过程

解：

原始数据序列： ${x_t}$ ，一阶差分变换后的序列为： $y_t$ 。

（1）通过下文中程序的运行，我们可以从运行结果中得到 $$ y_t = 1.253y_{t-1} - 0.3522y_{t-2} + \varepsilon_{t} + 0.5022 \varepsilon_{t-1} $$ （2）根据下文的运算结果得到未来10年的预测值分别为:： $$ 6419.44740352031 \ 6668.77039934881\ 6861.19145947359\ 7014.42501092823\ 7138.60914121919\ 7240.20495400936\ 7323.73568451267\ 7392.59199500092\ 7449.42827658915\ 7496.37548147182 $$

3.程序

求解的MATLAB程序如下：

clc, clear
format long g

% 定义列向量 xt，其中包含了原始数据的时间序列。
xt = [119, 142, 144, 150, 165, 168, 200, ...
    216, 218, 185, 173, 181, 208, 240, 254, ...
    235, 222, 243, 275, 288, 292, 309, 310, 327, ...
    316, 339, 379, 417, 460, 489, 525, 580, 682, ...
    853, 956, 1104, 1355, 1512, 1634, 1879, ...
    2287, 2939, 3923, 4854, 5576, 6079]';

% 进行一阶差分变换，即计算 xt 中相邻元素之间的差值，生成一个新的列向量 yt
% 一阶差分变换可以将非平稳时间序列转换为平稳时间序列。
yt = diff(xt);

% 拟合arma模型
m = armax(yt, [2, 1])
% armax函数接受两个参数：时间序列数据和模型阶数[p,q]，其中p是自回归项的数量，q是滞后误差项的数量

% 计算yt的10期预测值
ythat = forecast(m, yt, 10);

% 计算原始数据的10期预测值
% 这一行代码用于计算原始数据 xt 的预测值
% xt(end)表示原始数据的最后一个观测值，cumsum(ythat)表示将ythat中每个元素累加得到的新的列向量
xthat = xt(end) + cumsum(ythat)

4.结果

（1）我们可以从运行结果中得到 $$ y_t = 1.253y_{t-1} - 0.3522y_{t-2} + \varepsilon_{t} + 0.5022 \varepsilon_{t-1} $$ （2）未来10年的预测值分别为:： $$ 6419.44740352031 \ 6668.77039934881\ 6861.19145947359\ 7014.42501092823\ 7138.60914121919\ 7240.20495400936\ 7323.73568451267\ 7392.59199500092\ 7449.42827658915\ 7496.37548147182 $$

习题8.5

1. 题目要求

2.解题过程

解：

（1）序列时序图

记原始序列为 {${x_t}$} ,序列时序图如下图所示，时序图显示该序列大致具有12个周期变化，周期的长度为9或10年，下面使用周期 T=10行计算。

（2）差分平稳

对原序列做10步差分，消除季节趋势，得到序列 {$y_t$} ，其中， $y_t = x_{t+10}-x_t$，差分后序列图如下图所示。时序图显示差分后序列基本平稳了。

（3）模型拟合

根据差分后序列的自相关和偏自相关的性质，尝试拟合ARMA模型，拟合的ARMA (1,10) 模型较理想，并且通过了白噪声检验，说明低阶的ARMA模型不适合拟合这个序列。

（4）

求预测值。求得下两个年度的预测值为4310和3674。

3.程序

求解的MATLAB程序如下：

clc, clear
format long g

a = [269, 321, 585, 871, 1475, 2821, 3928, 5943, 4950, ...
    2577, 523, 98, 184, 279, 409, 2285, 2685, 3409, 1824, ...
    409, 151, 45, 68, 213, 546, 1033, 2129, 2536, 957, ...
    361, 377, 225, 360, 731, 1638, 2725, 2871, 2119, 684, ...
    299, 236, 245, 552, 1623, 3311, 6721, 4254, 687, 255, ...
    473, 358, 784, 1594, 1676, 2251, 1426, 756, 299, 201, ...
    229, 469, 736, 2042, 2811, 4431, 2511, 389, 73, 39, 49, ...
    59, 188, 377, 1292, 4031, 3495, 537, 105, 153, 387, 758, ...
    1307, 3465, 6991, 6313, 3794, 1836, 345, 382, 808, ...
    1388, 2713, 3800, 309, 2985, 3790, 674, 71, 80, 108, ...
    229, 399, 1132, 2432, 3575, 2935, 1537, 529, 485, 662, ...
    1000, 1520, 2657, 3396]';

n = length(a);
% 用MATLAB的plot函数绘制a的图像
plot(a, '.-')

% 使用for循环遍历从第11个到最后一个数据元素，并对前10个数据元素和当前数据元素进行差分计算得到一个新的列向量b
for i = 11:n
    b(i-10) = a(i) - a(i-10); % 进行季节差分变换
end
b = b';
figure, plot(b, '.-')

% 计算b的自相关性和偏自相关性
figure, subplot(121), autocorr(b)
subplot(122), parcorr(b)
% 对b序列进行模型拟合
cs = armax(b, [1, 10]); % 拟合模型
figure, myres = resid(cs, b); % 计算残差向量并画出残差的自相关函数图
% 拟合模型的残差向量myres
[h, p, st] = lbqtest(myres.outputdata, 'lags', [6, 12, 18]); % 进行LBQ检验
% 注意，上面outputdata一定要加上，不然会报错

bhat = forecast(cs, b, 2); % 计算b的2期预测值
ahat = a(end-9:end-8) + bhat % 求原始序列的2期预测值

4.结果

后两个年度的预测值为4310和3674

习题8.6

1. 题目要求

2.解题过程

解：

（1）对所给时间序列建模

首先对此序列进行观察分析。下图为数据曲线图：

可以看出具有指数上升趋势，因此，对确定性部分先拟合一个指数增长模型，即 $$ X_t = \mu_t + Y_t, \mu_t = R_1e^{r_1t} $$ 这里各季节依次编号为$t = 1,2,\dots,100$。

然后确定性趋势的拟合。为了能用线性最小二乘法估计参数$R_1$和$r_1$, $\mu_t = R_1 e^{r_1t}$两边取对数，得： $$ \ln \mu_t = \ln R_1 + r_1t $$ 利用观测数据求得$\hat{R}_1 = 12.6385,\hat{r}_1$ = 0.0162。剩余平方和为1683.5371。剩余序列$Y_t$如下图所示，可以认为是平稳的：

对剩余序列拟合ARMA模型。$Y_t$自相关与偏自相关如下图所示：

可初步断定$Y_t$的适应模型为AR模型，逐增加AR模型阶数进行拟合，其残差方差图如下图所示：

因此，合适的模型为AR (2)，即 $$ Y_t = \varphi_1Y_{t-1}+\varphi_2Y_{t-2} + a_t $$ 参数估计为$\hat{\varphi}_1 = 0.5451,\hat{\varphi}_2 = 0.2478$。

建立组合模型。最后要以已估计出来的$ R_1,r_1,\varphi_1,\varphi_2 $的值为初始值用非线性最小二乘法对模型参数进行整体估计，模型整体可写为 $$ X_t = \mu_t + Y_t = R_1e^{r_1t} + \varphi_1(X_{t-1}-R_1e^{r_1(t-1)} )+\varphi_2(X_{t-2})-R_1e^{r_1(t-2)}+a_t $$ 最终的参数整体估计为 $$ \hat{R}_1=12.1089,\hat{r}_1=0.017,\hat{\varphi}_1 = 0.517, \hat{\varphi}_2 = 0.2397 $$ 残差平方和为739.4402。

（2）

对所给的时间序列进行两年（8个季度）的预报。用所建的模型以1970年第4几度即t = 100为原点进行预测，结果如下表所示。

l	0	1	2	3	4	5	6	7	8
t+l	t	t+1	t+2	t+3	t+4	t+5	t+6	t+7	t+8
$\hat{X}_t(l)$	62.1	65.8298	66.8384	68.562	70.0083	71.4879	72.9238	74.3507	75.768

3.程序

（1）

求解的MATLAB程序如下：

clc, clear

% 将数据按照每年的每个季度依次写入
a = [7.5, 8.9, 11.1, 13.4, 15.5, 15.7, 15.6, 16.7, 18, 17.4, 17.9, ...
    18.8, 17.6, 17, 16.1, 15.7, 15.9, 17.9, 20.3, 20.4, 20.2, 20.5, ...
    20.9, 20.9, 21.1, 21.4, 18.2, 20.1, 21.4, 21.3, 21.9, 21.3, ...
    20.4, 20.4, 20.7, 20.7, 20.9, 23, 24.9, 26.5, 25.6, 26.1, 27, ...
    27.2, 28.1, 28, 29.1, 28.3, 25.7, 24.5, 24.4, 25.5, 27, 28.7, ...
    29.1, 29, 29.6, 31.2, 30.6, 29.8, 27.6, 27.7, 29, 30.3, 31, 32.1, ...
    33.5, 33.2, 33.2, 33.8, 35.5, 36.8, 37.9, 39, 41, 41.6, 43.7, ...
    44.4, 46.6, 48.3, 50.2, 52.1, 54, 56, 53.9, 55.6, 55.4, 56.2, ...
    57.9, 57.3, 58.8, 60.4, 63.1, 83.5, 64.8, 65.7, 64.8, 65.6, 67.2, 62.1]';

n = length(a);
t0 = [46:1 / 4:71 - 1 / 4];
t = [1:100]';
xishu = [ones(n, 1), t];
cs = xishu \ log(a)
cs(1) = exp(cs(1))
ahat = cs(1) * exp(cs(2)*t);
cha = a - ahat
res = sum(cha.^2)
subplot(121), plot(t0, a, '*-')
subplot(122), plot(t0, cha, '.-')
figure, subplot(121), autocorr(cha)
subplot(122), parcorr(cha)
for i = 1:10
    cs2 = ar(cha, i);
    cha2 = resid(cs2, cha);
    myvar(i) = sum(cha2.outputdata.^2) / (100 - i);
end
figure, plot(myvar, '*-')

（2）

求解的MATLAB程序如下：

clc, clear

% 定义一个函数句柄 xt，它的输入参数是一个向量 cs 和一个矩阵 x。 
% x 矩阵的第一列是 a 向量的第二个元素到倒数第二个元素，第二列是 a 向量的第一个元素到倒数第三个元素，
% 第三列是一个列向量，它包含数字3到100
% 这些数字将用于预测未来的季度。
% 函数的输出是一个向量，表示用于预测季度的预测值。
xt = @(cs, x)cs(1) * (exp(cs(2)*x(:, 3)) - cs(3) * exp(cs(2)*(x(:, 3) - 1)) - ...
    cs(4) * exp(cs(2)*(x(:, 3) - 2))) + cs(3) * x(:, 1) + cs(4) * x(:, 2);
cs0 = [12.6385, 0.0162, 0.5451, 0.2478]';

% 将数据按照每年的每个季度依次写入
a = [7.5, 8.9, 11.1, 13.4, 15.5, 15.7, 15.6, 16.7, 18, 17.4, 17.9, ...
    18.8, 17.6, 17, 16.1, 15.7, 15.9, 17.9, 20.3, 20.4, 20.2, 20.5, ...
    20.9, 20.9, 21.1, 21.4, 18.2, 20.1, 21.4, 21.3, 21.9, 21.3, ...
    20.4, 20.4, 20.7, 20.7, 20.9, 23, 24.9, 26.5, 25.6, 26.1, 27, ...
    27.2, 28.1, 28, 29.1, 28.3, 25.7, 24.5, 24.4, 25.5, 27, 28.7, ...
    29.1, 29, 29.6, 31.2, 30.6, 29.8, 27.6, 27.7, 29, 30.3, 31, 32.1, ...
    33.5, 33.2, 33.2, 33.8, 35.5, 36.8, 37.9, 39, 41, 41.6, 43.7, ...
    44.4, 46.6, 48.3, 50.2, 52.1, 54, 56, 53.9, 55.6, 55.4, 56.2, ...
    57.9, 57.3, 58.8, 60.4, 63.1, 83.5, 64.8, 65.7, 64.8, 65.6, 67.2, 62.1]';

% 创建一个矩阵 x，包含3列，用于作为函数 xt 的输入参数
% 第一列是 a 向量的第二个元素到倒数第二个元素，第二列是 a 向量的第一个元素到倒数第三个元素，
% 第三列是一个列向量，它包含数字3到100，这些数字将用于预测未来的季度
x = [a(2:end-1), a(1:end-2), [3:100]'];
cs = lsqcurvefit(xt, cs0, x, a(3:end))
res = a(3:end) - xt(cs, x);
Q = sum(res.^2)
autocorr(res)
xhat = a;

for j = 101:108
    xhat(j) = cs(1) * (exp(cs(2)*j) - cs(3) * exp(cs(2)*(j - 1)) - ...
        cs(4) * exp(cs(2)*(j - 2))) + cs(3) * xhat(j-1) + cs(4) * xhat(j-2);
end

xhat101_108 = xhat(101:108)

4.结果

（1）结果见上文解题过程

（2）

l	0	1	2	3	4	5	6	7	8
t+l	t	t+1	t+2	t+3	t+4	t+5	t+6	t+7	t+8
$\hat{X}_t(l)$	62.1	65.8298	66.8384	68.562	70.0083	71.4879	72.9238	74.3507	75.768

l	0	1	2	3	4	5	6	7	8
t+l	t	t+1	t+2	t+3	t+4	t+5	t+6	t+7	t+8
\(\hat{X}_t(l)\)	62.1	65.8298	66.8384	68.562	70.0083	71.4879	72.9238	74.3507	75.768