时间序列分析

时间序列分析是一种数据分析方法，它研究的对象是代表某一现象的一串随时间变化而又相关联的数据系列，从而描述和探索该现象随时间发展变化的规律性。
时间序列分析方法包括：
(1)确定性时序分析
(2)随机性时序分析

平稳时序分析主要通过建立:
自回归模型( ，Autoregressive Models)、
滑动平均模型( ，Moving Average Models)
自回归滑动平均模型( ，Autoregressive Moving Average Models)
分析平稳的时间序列的规律，一般的分析程序可用下面框图表示：

1.自回归模型AR§
如果时间序列Xt是它的前期值和随机项的线性函数，即可表示为

称为P阶自回归模型，记为AR（P）
注1：实参数 称为自回归系数，是待估参数。随机项 at是相互独立的白噪声序列，且服从均值为0、方差为 的正态分布。随机项与滞后变量不相关。
注2：一般假定Xt均值为0，否则令

2．滑动平均模型 MA（q)
如果时间序列Xt是它的当期和前期的随机误差项的线性函数，即可表示为

称为q阶移动平均模型，记为MA（q）注：实参数 为移动平均系数，是待估参数。

引入滞后算子，并令
则模型可简写为

注1：移动平均过程无条件平稳。
注2：滞后多项式 的根都在单位圆外时，AR过程与MA过程能相互表出，即过程可逆，
,即为MA过程的逆转形式，也就是MA过程等价于无穷阶的AR过程。
注3：
满足平稳条件时， AR过程等价于无穷阶的MA 过程，即
3. 自回归移动平均模型ARMA(p,q)
如果时间序列Xt是它的当期和前期的随机误差项以及前期值的线性函数，即可表示为

称为（p,q）阶的自回归移动平均模型，记为ARMA(p,q)

在实际的社会经济现象中收集到的时序大多数是呈现出明显的趋势性或周期性，要把趋势和波动综合考虑进来，于是Xt = μt+Yt。其中μt 表示Xt中随时间变化的均值(往往是趋势值)，Yt是Xt中剔除μt后的剩余部分，表示零均值平稳过程。

Matlab有关时间序列分析命令 (matlab 2013)

在工具箱中可以通过命令garchset指定模型的结构，garchset的语法为
Spec=garchset(‘属性1’，属性1的值，’属性2’，属性2的值，…)
在工具箱中可以通过命令garchfit对模型中的参数进行估计，garchfit的语法为
[Coeff,Errors,LLF,Innovations,Sigmas]=garchfit(Spec,Series)
其中输入参数Spec指定模型的结构，Series 为时间序列的样本观测值。输出参数Coeff是模型的参数估计值，Errors是模型参数的标准差，LLF是最大似然估计法中的对数目标函数值，Innovations是残差向量，Sigmas是对应于Innovations 的标准差。

案例分析：

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

1、问题分析与模型建立

首先画出数据的走势图，这一时间序列是具有明显趋势且不含有周期性变化经济波动序列，即为非平稳的时间序列，对此序列进行建模预测需要用上面介绍的非平稳时间序列分析方法。采用模型：Xt = μt+Yt
2. 模型求解
2.1 确定性趋势
确定趋势是按指数趋势发展的

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] ;

a = exp(-290.4864)=6.9721×10-127
b = exp(0.1510)=1.1630
2.2 随机性趋势
残差序列

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

2.3ωt的时间序列分析
ωt为上述二次差分序列的均值化序列
%序列{ωt}的样本自相关函数ρk程序
w=r2-mean(r2); %零均值化
gamao=var(w); %求方差
for j=1:27
gama(j)=w(j+1:end)*w(1:end-j)’/27;
end
rho=gama/gamao %样本自相关系数
bar(rho) %条状图

%样本偏相关函数Φ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) %条形图


模型定阶程序：
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=547.624381,BIC=551.511892
R=0,M=2,AIC=547.846439,BIC=553.029786
R=0,M=3,AIC=549.327469,BIC=555.806653
R=1,M=0,AIC=551.953049,BIC=555.840559
R=1,M=1,AIC=550.238380,BIC=555.421727
R=1,M=2,AIC=821.672996,BIC=828.152180
R=1,M=3,AIC=547.002577,BIC=554.777598
R=2,M=0,AIC=555.026386,BIC=560.209733
R=2,M=1,AIC=556.460271,BIC=562.939455
R=2,M=2,AIC=552.958726,BIC=560.733748
R=2,M=3,AIC=552.882696,BIC=561.953554
R=3,M=0,AIC=554.379959,BIC=560.859143
R=3,M=1,AIC=553.753258,BIC=561.528280
R=3,M=2,AIC=560.830281,BIC=569.901139
R=3,M=3,AIC=829.833618,BIC=840.200313
计算结果显示可以为ARMA（1,3）模型。
对模型 进行参数估计程序：
spec = garchset(‘R’,1,‘M’,3,‘Display’,‘off’); %指定模型的结构
[coeffX,errorsX,LLFX] = garchfit(spec,w) %拟合参数

模型的检验和预测程序：
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说明模型是可用的，未来三年ωt的预测值为15270，-25983，24497。
通过指数预测的值加上残差的预测值并可以得到最终预测结果。