第八章：时间序列

最新推荐文章于 2023-10-18 00:52:48 发布

victor_cs_

最新推荐文章于 2023-10-18 00:52:48 发布

阅读量2.5k

点赞数 3

分类专栏：数学建模算法与应用文章标签：时间序列数学建模

本文链接：https://blog.csdn.net/qq_41000485/article/details/99592690

版权

数学建模算法与应用专栏收录该内容

16 篇文章 25 订阅

订阅专栏

8.1 确定型时间序列分析方法

8.1.1 移动平均法

只适合做近期预测，而且是目标的发展趋势变化不大的情况。

clc,clear
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,5];   %n为移动平均的项数
for i=1:length(n)    %由于n的取值不同，因此下面是用来细胞数组
    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)); %求预测的标准误差
end
y12, s   %分别显示两种方法的预测值和预测的标准误差

8.1.2 指数平滑法

（1）一次指数平滑法
加权系数 $\alpha$ 的选择：

时间序列波动不大，比较平稳，则 $\alpha$ 应取小一点
时间序列具有迅速且明显的变动趋向，则 $\alpha$ 应取大一点

clc,clear
yt=load('dianqi.txt');   %实际销售额数据以列向量的方式存放在纯文本文件中
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
yhat
err=sqrt(mean((repmat(yt,1,m)-yhat).^2)) 
xlswrite('dianqi.xls',yhat) %把预测数据写到Excel文件，准备在Word表格中使用
yhat1988=alpha*yt(n)+(1-alpha).*yhat(n,:)

（2）二次指数平滑法

clc,clear
yt=load('fadian.txt'); 
n=length(yt), alpha=0.3; st1(1)=yt(1); st2(1)=yt(1);
for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
end
xlswrite('fadian.xls',[st1',st2']) %把数据写入表达Sheet1中的前两列
at=2*st1-st2;
bt=alpha/(1-alpha)*(st1-st2);
yhat=at+bt;  %最后一个分量为1986年的预测值
xlswrite('fadian.xls',yhat','Sheet1','C2') %把预测值写入第三列
str=['C',int2str(n+2)]; %准备写1987年的预测值位置的字符串
xlswrite('fadian.xls',at(n)+2*bt(n),'Sheet1',str)%把1987年的预测值写到相应位置

（3）三次指数平滑法

clc,clear
yt=load('touzi.txt'); 
n=length(yt); alpha=0.3; st0=mean(yt(1:3));
st1(1)=alpha*yt(1)+(1-alpha)*st0;
st2(1)=alpha*st1(1)+(1-alpha)*st0;
st3(1)=alpha*st2(1)+(1-alpha)*st0;
for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
    st3(i)=alpha*st2(i)+(1-alpha)*st3(i-1);
end
xlswrite('touzi.xls',[st1',st2',st3'])  
at=3*st1-3*st2+st3;
bt=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1-2*(5-4*alpha)*st2+(4-3*alpha)*st3);
ct=0.5*alpha^2/(1-alpha)^2*(st1-2*st2+st3);
yhat=at+bt+ct;
xlswrite('touzi.xls',yhat','Sheet1','D2') 
plot(1:n,yt,'D',2:n,yhat(1:end-1),'*')
legend('实际值','预测值','Location','northwest')  %图注显示在左上角
xishu=[ct(end),bt(end),at(end)]; %二次预测多项式的系数向量
yhat1990=polyval(xishu,2)  %求预测多项式m=2时的值

8.1.3 差分指数平滑法

（1）一节差分指数平滑法
时间序列呈直线增加

clc,clear
yt=load('ranliao.txt'); 
n=length(yt); alpha=0.4; 
dyt=diff(yt); %求yt的一阶向前差分
dyt=[0;dyt]; %这里使用的是一阶向后差分，加0补位
dyhat(2)=dyt(2); %指数平滑值的初始值
for i=2:n
    dyhat(i+1)=alpha*dyt(i)+(1-alpha)*dyhat(i);
end
for i=1:n
    yhat(i+1)=dyhat(i+1)+yt(i);  
end
yhat
xlswrite('ranliao.xls',[yt,dyt]) 
xlswrite('ranliao.xls',[dyhat',yhat'],'Sheet1','C1')

（2）二阶差分指数平滑型
所产生的新序列基本上是平稳的，只能逐期预测

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

季节系数法

clc, clear
format long g
a=load('jijie.txt');
[m,n]=size(a);
a_mean=mean(mean(a));  %计算所有数据的算术平均值
aj_mean=mean(a);  %计算同季节的算术平均值
bj=aj_mean/a_mean  %计算季节系数
w=1:m;
yhat=w*sum(a,2)/sum(w)  %预测下一年的年加权平均值，这里是求行和
yjmean=yhat/n  %计算预测年份的季节平均值
yjhat=yjmean*bj  %预测年份的季节预测值
format   %恢复默认的显示格式

8.2 平稳时间序列模型

是指宽平稳，即均值和协方差不随时间的平移而变化

8.2.1 时间序列的基本概念

随机过程的均值函数
随机过程的方差函数
自相关函数是标准化自协方差函数
平稳序列
平稳白噪声序列
随机线性序列
偏相关函数

平稳时间序列——ARMA时间序列

AR（p）序列：自回归序列
MA（q）序列：移动平均序列
ARMA（p，q）序列：自回归移动平均序列

8.2.2 ARMA模型的构建及预报

AIC定阶准则
参数估计（使用Matlab工具箱）
卡方检验（使用lbqtest函数）

8.3 时间序列的Matlab相关工具箱及命令

8.3.1 系统辨识工具箱

观察数据：idinput
数据预处理：detrend，idfilt，idresamp
模型结构选择：struc，arxstruc，ivstruc，selstruc
参数估计：ar，arx，armax，ivx
模型预报与仿真：compare，pe，predict，sim

clc, clear
elps=randn(10000,1); x(1:2)=0;
for i=3:10000
    x(i)=-0.6*x(i-1)-0.2*x(i-2)+elps(i); %产生模拟数据
end
xlswrite('data1.xls',x(end-9:end)) 
dlmwrite('mydata.txt',x)  
x=x'; m=ar(x,2)   %进行参数估计
xhat=forecast(m,x,3) %计算3个预测值

clc, clear
elps=randn(10000,1); x(1:2)=0;
for i=3:10000
    x(i)=-0.6*x(i-1)-0.2*x(i-2)+elps(i);
end
for i=1:3
    m{i}=ar(x,i);   %拟合模型
    myaic(i)=aic(m{i}); %计算AIC的值
end
myaic

clc, clear
elps=randn(10000,1); x(1:2)=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);  %拟合模型
xp=predict(m,x);  %计算已知的数据值
res=x-xp;  %计算残差向量，也可以使用命令res=resid(m,x)计算残差向量
h=lbqtest(res) %对残差向量进行Ljung-Box检验

clc, clear
elps=randn(10000,1);
for i=3:10000
    x(i)=elps(i)-0.6*elps(i-1)-0.2*elps(i-2); 
end
x=x'; m=armax(x,[0,2])  %拟合ARMA(0,2)模型，x必须为列向量

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
p=input('输入阶数p=');q=input('输入阶数q='); %输入模型的阶数
m=armax(x,[p,q])  %拟合指定参数p，q的模型
res=resid(m,x); %计算残差向量
h=lbqtest(res) %进行chi2检验

8.3.2 计量经济学工具箱

模型参数获取或设置：garchget，garchset
建模和仿真：garchfit，garchpred，garchsim

CARCH模型——广义自回归条件异方差

clc,clear
VarMd=garch('Constant',0.01,'GARCH',0.2,'ARCH',0.3); %指定模型的结构
Md=arima('Constant',0,'AR',0.8,'MA',0.4, 'Variance',VarMd); %指定模型的结构
[y,e,v]=simulate(Md,10000); %产生指定结构模型10000个模拟数据
ToEstVarMd=garch(1,1);
ToEstMd=arima('ARLags',1,'MALags',1,'Constant',0,'Variance',ToEstVarMd);
[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,y)  %模型拟合
res=infer(EstMd,y); %计算残差
h=lbqtest(res) %进行模型检验
yhat=forecast(EstMd,3,'Y0',y) %预测未来的三个值

clc, clear
ToEstVarMd=garch(1,1);
%ToEstMd=arima('ARLags',1:2,'Variance',ToEstVarMd); %AR的阶次取2，无法通过
ToEstMd=arima('ARLags',1,'Variance',ToEstVarMd);
y=load('mydata.txt');
[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,y')  %模型拟合，注意y为列向量
yhat=forecast(EstMd,3,'Y0',y')

8.3.3 金融工具箱

构造时间序列数组：fints，ascii2fts
用户图形界面：ftsgui，ftstool

price=[1:6]'
dates=[today:today+5]'
obj=fints(dates,price)

8.4 ARIMA序列和季节性序列

8.4.1 ARIMA序列及其预报

clc,clear
a=textread('hua.txt');  %把原始数据按照原来的排列格式存放在纯文本文件中
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)); %计算拟合参数的个数
        %compute Akaike and Bayesian Information Criteria
        [aic(k),bic(k)]=aicbic(logL,numParams,n);
    end
end
fprintf('R,M,AIC,BIC的对应值如下\n %f');  %显示计算的结果
check=[R',M',aic',bic']
r=input('输入阶数p=');m=input('输入阶数q=');
ToEstMd=arima('ARLags',1:r,'MALags',1:m,'Constant',0); %指定模型的结构
[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,da); %模型拟合
dx_Forecast=forecast(EstMd,10,'Y0',da)  %计算10步预报值
x_Forecast=a(end)+cumsum(dx_Forecast)   %计算原始数据的10步预测值

取简化模型IMA(1，1)

clc,clear
a=textread('hua.txt');  
a=nonzeros(a'); 
da=diff(a);      %计算一阶差分
ToEstMd=arima('MALags',1); 
[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,da); 
dx_Forecast=forecast(EstMd,10,'Y0',da)  
x_Forecast=a(end)+cumsum(dx_Forecast)

8.4.2 季节性序列及其预报

clc,clear
x=load('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);   %消除趋势性的差分运算
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)); %计算拟合参数的个数
        %compute Akaike and Bayesian Information Criteria
        [aic(k),bic(k)]=aicbic(logL,numParams,m2);
    end
end
fprintf('R,M,AIC,BIC的对应值如下\n %f');  
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')  
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个预报值

victor_cs_

关注

3
点赞
踩
37

收藏

觉得还不错? 一键收藏
2
评论
第八章：时间序列

8.1 确定型时间序列分析方法8.1.1 移动平均法只适合做近期预测，而且是目标的发展趋势变化不大的情况。clc,cleary=[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,5]; %n为移动平均的项数for i=1:leng...
复制链接

扫一扫