导入每年的太阳黑子数数据loadsunspotdatwhosyearsunspot.docx
《导入每年的太阳黑子数数据loadsunspotdatwhosyearsunspot.docx》由会员分享，可在线阅读，更多相关《导入每年的太阳黑子数数据loadsunspotdatwhosyearsunspot.docx（13页珍藏版）》请在冰豆网上搜索。
导入每年的太阳黑子数数据loadsunspotdatwhosyearsunspot
1.导入每年的太阳黑子数数据
loadsunspot.dat
whos
year=sunspot（:
1）;
wolfer=sunspot（:
2）;
plot（year,wolfer）
xlabel（'Year'）
ylabel（'Numberofsunspot'）
title（'SunspotData'）
考察序列的图形，该序列为近似的平稳序列或含有周期大约为11的季节项。
2．计算样本自相关函数和偏相关函数
figure
autocorr（wolfer,40,2）;
legend（'自相关系数函数'）;
说明：
可认为该序列的自相关系数函数随时间振荡缓慢而趋向于0；
可认为该序列具有周期为11的特征，建议先消去周期性趋势。
Syntax
[PartialACF,Lags,Bounds]=parcorr（Series,nLags,R,nSTDs）
Syntax
[PartialACF,Lags,Bounds]=parcorr（Series,nLags,R,nSTDs）autocorr（Series,nLags,M,nSTDs）
[ACF,Lags,Bounds]=autocorr（Series,nLags,M,nSTDs）
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（'Numberofsunspot'）
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）;
forp=0:
3
forq=0:
3
m=armax（wolfer0,[p,q]）;
a（p+1,q+1）=aic（m）;
end
end
a=
6.99786.15955.68545.5724
5.88225.64345.53375.5249
5.47965.48195.48275.4920
5.48395.48365.40025.3914
mw=armax（wolfer0,[33]）;
mw
Discrete-timeIDPOLYmodel:
A（q）y（t）=C（q）e（t）
A（q）=1-2.497q^-1+2.338q^-2-0.8131q^-3
C（q）=1-1.427q^-1+0.5083q^-2+0.09392q^-3
EstimatedusingARMAXfromdatasetwolfer0
Lossfunction210.311andFPE219.521
Samplinginterval:
1
模型方程：
5.所建模型的白噪声检验
（1）残差计算及正态分布检验法
>>e=resid（mw,wolfer0）;
图所表示的是残差的
数据图（The99%confidenceintervalsforthesevaluesarealsocomputedandshownasayellowregion），没有一个超出其范围，不能否定为白噪声。
（2）白噪声的
检验
rho=autocorr（e,40）;
form=1:
15
sum2=0;
fork=1:
m
sum2=sum2+rho（k+1）*rho（k+1）;
end
chipf（m）=288*sum2;
lamda（m）=chi2inv（0.95,m）;
end
t=1:
15;
plot（t,chipf（t）,t,lamda（t））;
legend（'chipf','临界值lamda'）;
（3）递推预测残差图
plot（e）;
6.模型递推预测值与实际观测值比较图
wolferpredict=predict（mw,wolfer0）;
plot（year,wolfer0）;
holdon
plot（year,wolferpredict,'r'）;
holdoff
xlabel（'Year'）
ylabel（'Numberofsunspot'）
title（'modelpredictnumberVersusExactnumber'）
注：
也可利用garchToolbox建模
spec=garchset（'R',2,'M',1,'P',0,'Q',0,'Display','off'）;
[spec,errors,LLF,residuals,sigmas]=garchfit（spec,wolfer0）;
garchdisp（spec,errors）%Displaytheestimationresults
plot（year（1:
end）,residuals）
xlabel（'Date'）,ylabel（'Residual'）,title（'FilteredResiduals'）
Mean:
ARMAX（2,1,0）;Variance:
GARCH（0,0）
ConditionalProbabilityDistribution:
Gaussian
NumberofModelParametersEstimated:
5
StandardT
ParameterValueErrorStatistic
---------------------------------------------
C0.0401220.959930.0418
AR
（1）1.37810.07423718.5629
AR
（2）-0.645490.067808-9.5195
MA
（1）-0.11550.094943-1.2165
K232.7415.7814.7494
standardizedResiduals=residuals./sigmas;
figure
autocorr（standardizedResiduals,25）
title（'SampleACFofStandardizedResiduals'）
figure
parcorr（standardizedResiduals,25）
title（'SamplePACFofStandardizedResiduals'）
升级会员 