1 导入每年的太阳黑子数数据 load sunspotdat whos year=sunspot.DOC

1. 导入每年的太阳黑子数数据 load sunspot.dat whos year=sunspot(:,1); wolfer=sunspot(:,2); plot(year,wolfer) xlabel(Year) ylabel(Number of sunspot) title(Sunspot Data) 考察序列的图形，该序列为近似的平稳序列或含有周期大约为11的季节项。 2．计算样本自相关函数和偏相关函数 figure autocorr(wolfer,40,2); legend(自相关系数函数); 说明：可认为该序列的自相关系数函数随时间振荡缓慢而趋向于0； 可认为该序列具有周期为11的特征，建议先消去周期性趋势。 figure parcorr(wolfer,40); legend(偏相关系数函数); 说明：可认为该序列的偏相关系数函数随时间振荡缓慢而趋向于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)); 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*ones(1,27); wolfer3=reshape(wolfer2,1,297); wolfer0=wolfer3(1:288); plot(year,wolfer0); （3）计算样本自相关函数和偏相关函数 figure autocorr(wolfer0,40,2); legend(消去季节项后的数据的自相关系数函数); figure parcorr(wolfer0,40); legend(消去季节项后的数据的偏相关系数函数); 样本自相关函数和偏相关函数的特征与原来变化不大，随时间振荡趋向于0，不具有较明显的截尾性质，建议采ARMA模型进行建模。 4. ARMA(p,q)模型. 对 p=0,1,2,3； q=0,1,2,3; 逐对建立模型，取AIC达到最小的为合适模型 a=zeros(4,4); for p=0:3 for q=0:3 m=armax(wolfer0,[p,q]); a(p+1,q+1)=aic(m); end end a = 6.9978 6.1595 5.6854 5.5724 5.8822 5.6434 5.5337 5.5249 5.4796 5.4819 5.4827 5.4920 5.4839 5.4836 5.4002 5.3914 mw= armax(wolfer0,[3 3]); mw Discrete-time IDPOLY model: A(q)y(t) = C(q)e(t) A(q) = 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