matlab实习实验报告

实验一时间序列ARMA学习

一、基本形式

ARMA模型分为以下三种：

1.自回归模型（AR）

如果时间序列满足

其中是独立同分布的随机变量序列，且满足：

以及E( ) = 0，则称时间序列为服从p阶的自回归模型。

AR模型是用过去的观测值和现在的干扰值的线性组合来预测未来值，式子中的p为阶数，β为待定参数，y为一个平稳序列。

自回归模型的平稳条件：

滞后算子多项式

的根均在单位圆外，即φ(B) = 0的根大于1。

2.移动平均模型（MA）

如果时间序列满足

则称时间序列为服从q阶移动平均模型；

MA模型是用过去的干扰值和现在的干扰值的线性组合预测未来的值，式中q为阶数，α为待定参数，y为一个平稳序列。

移动平均模型平稳条件：任何条件下都平稳。

3.自回归滑动平均模型（ARMA）

如果时间序列满足：

则称时间序列为服从(p,q)阶自回归滑动平均混合模型。

模型求解流程图：

二、实例分析

太阳的光球表面有时会出现一些暗的区域，它是磁场聚集的地方，这就是太阳黑子。黑子是太阳表面可以看到的最突出的现象。一个中等大小的黑子大概和地球的大小差不多。长期的观测发现，黑子多的时候，其他太阳活动现象也会比较频繁。从外文网站查询到太阳黑子数量的变化，从而预测接下来数量的走势。

首先整理数据，进行平稳性判断，使用matlab语言中的adftest()函数即可。

数据为平稳数据，因而无需做平稳化处理，使用AR模型对原数据作出预测。

Matlab代码为：

clc

clear

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)%

figure(2);

subplot(2,1,1)

autocorr(c,500)%²

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50);

set(gca,'Xlim',[0 50]);

z=iddata(c);%×

m=armax(z(1:end),'na',20,'nc',0);%

yp = predict(m,c,1);

figure(3)

plot(c,'-.');

hold on

plot(yp,'r')

grid

legend('³õʼֵ','Ô¤²âÖµAR')

bzc=sum((yp-c).^2)

结果如下：

2.MA模型求解

Matlab代码：

clc

clear

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)

figure(2);

subplot(2,1,1)

autocorr(c,500)

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50);

set(gca,'Xlim',[0 50]);

z=iddata(c);

m=armax(z(1:end),'na',0,'nc',4); yp = predict(m,c,1);

figure(3)

plot(c,'-.');

hold on

plot(yp,'r')

grid

legend('³õʼֵ','Ô¤²âÖµMA')

bzc=sum((yp-c).^2)

结果如下：

3.ARMA模型求解

对于ARMA模型，需要做出阶数判断，用ACF和PACF确定p、q阶数,在使用ARMA模型求解。

Matlab代码：

clc

clear

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)

figure(2);

subplot(2,1,1)

autocorr(c,500)%²自相关

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50); %偏自相关

set(gca,'Xlim',[0 50]);

结果如下：

由图可以看出ACF和PACF都是拖尾，故可用ARMA模型求解，模拟数据matlab代码如下：

a=csvread('heizi.csv');

b=a(2700:3234,[1,4]);

c=b(:,2);

diff(c);

adftest(c)%

figure(2);

subplot(2,1,1)

autocorr(c,500)%

set(gca,'Xlim',[0 500]);

subplot(2,1,2)

parcorr(c,50); %

set(gca,'Xlim',[0 50]);

z=iddata(c);

m=armax(z(1:end),'na',19,'nc',19);

yp = predict(m,c,1);

figure(3)

plot(c,'-.');

hold on

plot(yp,'r')