用时间序列分析的知识解决国内生产总值序列问题.docx

国内生产总值序列分析

一．问题的提出
选取1978-2006历年国内生产总值数据如下，试对该时间序列进行建模并预测。

二．问题分析与模型建立
首先画出数据的走势图，这一时间序列是具有明显趋势且不含有周期性变化经济波动序列，即为非平稳的时间序列，对此序列进行建模预测需要用上面介绍的非平稳时间序列分析方法。采用模型：
Xt=μt+Yt
其中μt表示Xt中随时间变化的趋势值，Yt是Xt中剔除μt后剩余部分。

历年国内生产总值时间序列图

三．模型求解
1．确定性趋势
确定趋势是按指数趋势发展的μt=abt------àlnμt=lna+tlnb
线性回归分析程序：
t=1978:2006;
x=[3624.10 4038.20 4517.80 4862.40 5294.70 5934.50 7171.00 8964.40 10202.20 11962.50 14928.30 16909.20 18547.90 21617.80 26638.10 34634.40 46759.40 58478.10 67884.60 74462.60 78345.20 82067.46 89468.10 97314.80 105172.34 116898.40 136515.00 182321.00 209407.00];
X=[ones(29,1) t']; %回归的资料矩阵
y=log(x)'; %线性化
[B,BINT,R,RINT,STATS] = regress(y,X) %回归
y2= exp(B(1)+B(2).*t) %预测值
plot(t,x,t,y2,'+'); %回归效果图

>>B =-290.4864
0.1510
STATS = 1.0e+003 *
0.0010 2.1838 0 0.0000

原始数据与指数回归数据对比图
得到B=[-290.4864,0.1510]
STATS=1.0e+003*[0.0010,2.1838,0]
即

由上图可知仅用指数回归的效果较差。
2.随机性趋势
（1）残差序列Yt={Xt-μt}
r=x-y2; %残差数列
plot(t,r,'O'); %残差散点图

参差序列散点图
观察残差序列的散点图可知，该序列有很大的波动性，可认为是非平稳的，应该经过多次差分使其平稳。
（2）二次差分后序列∇Yt=Yt-2Yt-1+Yt-2
r1=diff(r); %残差的一阶差分
r11=[0 r1]; %补数列差分后的项为0
plot(t,r11,'o'); %一阶差分散点图
r2=diff(r1); %二阶差分
r21=[0 0 r2]; %补数列差分后的项为0
plot(t,r21,'o'); %二阶差分散点图

一阶差分散点图

二阶差分散点图
（3）wt的时间序列分析
A. 将序列{r2t}零均值化,序列{wt}的样本自相关函数pk程序如下：
w=r2-mean(r2); %零均值化
gamao=varw); %求方差
for j=1:27
gama(j)=w(j+1:end)*w(1:end-j)'/27;
end
rho=gama/gamao %样本自相关系数
bar(rho) %条状图

自相关系数条形图
B. 样本偏相关函数φkk程序如下：
f(1,1)=rho(1);
for k=2:27
s1=rho(k);s2=1; %计算的初始值
for j=1:k-1
s1=s1-rho(k-j)*f(k-1,j);
s2=s2-rho(j)*f(k-1,j);
end
f(k,k)=s1/s2; %对角上的样本偏相关系数
for j=1:k-1
f(k,j)=f(k-1,j)-f(k,k)*f(k-1,k-j); %不在对角上的样本偏相关系数
end
end
pcorr=diag(f)' %提取偏相关函数
bar(pcorr) %条形图

偏自相关函数

C. 模型定阶的程序：
for i=0:3
for j=0:3
spec= garchset('R',i,'M',j,'Display','off'); %指定模型的结构
[coeffX,errorsX,LLFX] = garchfit(spec,w); %拟合参数
num=garchcount(coeffX); %计算拟合参数的个数
[aic,bic]=aicbic(LLFX,num,27);
fprintf('R=%d,M=%d,AIC=%f,BIC=%f\n',i,j,aic,bic); %显示计算结果
end
end
结果如下：
R=0,M=0,AIC=554.744695,BIC=557.336369
R=0,M=1,AIC=548.981658,BIC=552.869169
R=0,M=2,AIC=548.671841,BIC=553.855188
R=0,M=3,AIC=550.112192,BIC=556.591376
R=1,M=0,AIC=550.968125,BIC=554.855636
R=1,M=1,AIC=550.239945,BIC=555.423293
R=1,M=2,AIC=551.360349,BIC=557.839534
R=1,M=3,AIC=546.975261,BIC=554.750283
R=2,M=0,AIC=552.918590,BIC=558.101938
R=2,M=1,AIC=559.057147,BIC=565.536332
R=2,M=2,AIC=551.171163,BIC=558.946184
R=2,M=3,AIC=552.530372,BIC=561.601230
R=3,M=0,AIC=553.140182,BIC=559.619367
R=3,M=1,AIC=553.153087,BIC=560.928109
R=3,M=2,AIC=561.426269,BIC=570.497127
R=3,M=3,AIC=553.156375,BIC=563.523070
得到结果显示，可以认为是ARMA（1,3）模型
D. 对模型wt=c+φ1wt-1+εt+θ1εt-1+θ2εt-2+θ3εt-3
进行参数估计程序：
spec = garchset('R',1'M',3,'Display','off'); %指定模型的结构
[coeffX,errorsX,LLFX] = garchfit(spec,w) %拟合参数
运行结果如下:
coeffX =
Comment: 'Mean: ARMAX(1,3,0); Variance: GARCH(0,0)'
Distribution: 'Gaussian'
R: 1
M: 3
C: 35.8925
AR: -0.6250
MA: [-0.0387 0.0387 -1.0000]
VarianceModel: 'GARCH'
K: 3.5044e+007
Display: 'off'
于是ARMA（1,3）模型为
wt=35.8925-0.6250wt-1+εt-0.0387εt-1+0.0387εt-2-1.0000εt-3
E. 模型的检验和预测程序：
spec= garchset('R',1,'M',3); %指定模型的结构
[coeff,errors,LLF,innovations,sigmas,summary] = garchfit(spec,w) %拟合参数
h=lbqtest(innovations) %模型检验
[sigmaForecast,x_Forecast] = garchpred(coeff,w,3) %预测
得到的结果为 h=0，说明模型是可用的，未来三年 wt的预测值为 12010 -24258 20979
（4） 通过指数预测的值加上残差的预测值并可以得到最终预测结果。

7 / 7