导入每年的太阳黑子数数据loadsunspotdatwhosyearsunspot.docx

1、导入每年的太阳黑子数数据loadsunspotdatwhosyearsunspot1. 导入每年的太阳黑子数数据load sunspot.datwhosyear=sunspot(:,1);wolfer=sunspot(:,2);plot(year,wolfer)xlabel(Year)ylabel(Number of sunspot)title(Sunspot Data)考察序列的图形，该序列为近似的平稳序列或含有周期大约为11的季节项。2计算样本自相关函数和偏相关函数figureautocorr(wolfer,40,2);legend(自相关系数函数);说明：可认为该序列的自相关系数函数随时

2、间振荡缓慢而趋向于0； 可认为该序列具有周期为11的特征，建议先消去周期性趋势。SyntaxPartialACF,Lags,Bounds = parcorr(Series,nLags,R,nSTDs)SyntaxPartialACF,Lags,Bounds = parcorr(Series,nLags,R,nSTDs)autocorr(Series,nLags,M,nSTDs)ACF,Lags,Bounds = autocorr(Series,nLags,M,nS TDs)figureparcorr(wolfer,40);legend(偏相关系数函数);说明：可认为该序列的偏相关系数函数随时间

3、振荡缓慢而趋向于0 3. 时间序列的分解（消去季节项）（1）季节项的估计（取每个周期内的数据平均值为各自的季节项估计）；length(wolfer);wolfer1=wolfer; zeros(9,1);wolferm=reshape(wolfer1,11,27);meanm=mean(wolferm(1:2,1:27) mean(wolferm(3:11,1:26);（等价于meanm=mean(wolferm(1,1:27);mean(wolferm(2,1:27);mean(wolferm(3,1:26);mean(wolferm(4,1:26);mean(wolferm(5,1:26)

4、;mean(wolferm(6,1:26);mean(wolferm(7,1:26);mean(wolferm(8,1:26);mean(wolferm(9,1:26);mean(wolferm(10,1:26);mean(wolferm(11,1:26);）meanm1= meanm*ones(1,27);seasonal=meanm1(1:288);plot(year,seasonal);xlabel(Year)ylabel(Number of sunspot)title(太阳黑子数的季节项估计)（2）消去季节项后数据,其样本均值可近似认为为0；wolfer2=wolferm-meanm*

5、ones(1,27);wolfer3=reshape(wolfer2,1,297);wolfer0=wolfer3(1:288);plot(year,wolfer0);（3）计算样本自相关函数和偏相关函数figureautocorr(wolfer0,40,2);legend(消去季节项后的数据的自相关系数函数);figureparcorr(wolfer0,40);legend(消去季节项后的数据的偏相关系数函数);样本自相关函数和偏相关函数的特征与原来变化不大，随时间振荡趋向于0，不具有较明显的截尾性质，建议采ARMA模型进行建模。4. ARMA(p,q)模型. 对 p=0,1,2,3； q=

6、0,1,2,3; 逐对建立模型，取AIC达到最小的为合适模型a=zeros(4,4);for p=0:3for q=0:3m=armax(wolfer0,p,q);a(p+1,q+1)=aic(m);endenda = 6.9978 6.1595 5.6854 5.5724 5.8822 5.6434 5.5337 5.5249 5.4796 5.4819 5.4827 5.49205.4839 5.4836 5.4002 5.3914mw= armax(wolfer0,3 3);mwDiscrete-time IDPOLY model: A(q)y(t) = C(q)e(t) A(q) =

7、1 - 2.497 q-1 + 2.338 q-2 - 0.8131 q-3 C(q) = 1 - 1.427 q-1 + 0.5083 q-2 + 0.09392 q-3 Estimated using ARMAX from data set wolfer0 Loss function 210.311 and FPE 219.521 Sampling interval: 1 模型方程：5. 所建模型的白噪声检验（1）残差计算及正态分布检验法 e=resid(mw,wolfer0); 图所表示的是残差的数据图（The 99% confidence intervals for these val

8、ues are also computed and shown as a yellow region），没有一个超出其范围，不能否定为白噪声。（2）白噪声的检验rho=autocorr(e,40);for m=1:15sum2=0;for k=1:msum2=sum2+ rho(k+1)* rho(k+1);endchipf(m)=288*sum2;lamda(m)=chi2inv(0.95,m);endt=1:15;plot(t,chipf(t),t,lamda(t);legend(chipf,临界值lamda);（3）递推预测残差图plot(e);6. 模型递推预测值与实际观测值比较图wo

9、lferpredict = predict(mw,wolfer0);plot(year,wolfer0);hold onplot(year,wolferpredict,r);hold offxlabel(Year)ylabel(Number of sunspot)title(model predict number Versus Exact number)注：也可利用garch Toolbox建模spec = garchset(R,2,M,1 ,P, 0, Q, 0,Display,off); spec, errors, LLF, residuals, sigmas = garchfit(sp

10、ec, wolfer0);garchdisp(spec, errors) % Display the estimation resultsplot(year(1:end), residuals)xlabel(Date), ylabel(Residual), title (Filtered Residuals)Mean: ARMAX(2,1,0); Variance: GARCH(0,0) Conditional Probability Distribution: Gaussian Number of Model Parameters Estimated: 5 Standard T Parame

11、ter Value Error Statistic - - - - C 0.040122 0.95993 0.0418 AR(1) 1.3781 0.074237 18.5629 AR(2) -0.64549 0.067808 -9.5195 MA(1) -0.1155 0.094943 -1.2165 K 232.74 15.78 14.7494standardizedResiduals = residuals ./ sigmas;figureautocorr(standardizedResiduals,25)title(Sample ACF of Standardized Residuals)figureparcorr(standardizedResiduals,25)title(Sample PACF of Standardized Residuals)