某化工生产过程每2小时的浓度读数如表13所示。

表13  化工生产过程浓度数据

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.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

经计算破10步预报值见表14。

表14  10步预报值

步数 1 2 3 4 5

预报值 17.3611 17.4446 17.5811 17.6441 17.5821

步数 6 7 8 9 10

预报值 17.4644 17.4079 17.4644 17.5759 17.6364

clc,clear

a=textread('hua.txt');  %把原始数据按照原来的排列格式存放在纯文本文件hua.txt

a=nonzeros(a')'; %按照原来数据的顺序去掉零元素

r11=autocorr(a)   %计算自相关函数

r12=parcorr(a)   %计算偏相关函数

da=diff(a);      %计算1阶差分

r21=autocorr(da)  %计算自相关函数

r22=parcorr(da)   %计算偏相关函数

n=length(da);  %计算差分后的数据个数

for i=0:3

    for j=0:3

    spec= garchset('R',i,'M',j,'Display','off'); %指定模型的结构

    [coeffX,errorsX,LLFX] = garchfit(spec,da);  %拟合参数

    num=garchcount(coeffX);   %计算拟合参数的个数

    %compute Akaike and Bayesian Information Criteria

    [aic,bic]=aicbic(LLFX,num,n);  

    fprintf('R=%d,M=%d,AIC=%f,BIC=%f\n',i,j,aic,bic);  %显示计算结果

    end

end

r=input('输入阶数R');m=input('输入阶数M');

spec2= garchset('R',r,'M',m,'Display','off'); %指定模型的结构

[coeffX,errorsX,LLFX] = garchfit(spec2,da)    %拟合参数

[sigmaForecast,w_Forecast] = garchpred(coeffX,da,10)  %计算10步预报值

x_pred=a(end)+cumsum(w_Forecast)   %计算原始数据的10步预测值