- 积分
- 7031
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-1-1
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
AR预测模型
% Akaike信息检验准则率定模型阶数p
% 2013/6/14
% mingbo1501
%%
clear,clc,close all;
p=12;%模型阶数
num=500;%检验个数
%----------------生成样本--------------------%
data0=load('runoff.txt');
[m,n]=size(data0);
data=zeros(m-p,p+1);
for i=1:m-p
data(i,:)=data0(i:i+p);%前p列输入,第p+1列输出
end
[m1,n1]=size(data);
X1=data(1:m1-num,1:p);X1=[ones(m1-num,1),X1];
Y1=data(1:m1-num,p+1);%训练样本
X2=data(m1-num+1:end,1:p);X2=[ones(num,1),X2];
Y2=data(m1-num+1:end,p+1);%检验样本
%-----------------估计参数--------------------%
fai=(X1'*X1)\X1'*Y1;%最小二乘估计
yuce1=X1*fai;
err1=Y1-X1*fai;%残差
sigma=var(err1);%残差的方差
figure
plot(1:m1-num,yuce1,'b',1:m1-num,Y1,'r');
legend('预测值','实测值');
rate1=length(find((yuce1-Y1)./Y1*100<10))/(m1-num)*100;%合格率
DC1=1-sum((Y1-yuce1).^2)./sum((Y1-mean(Y1)).^2);%确定性系数
%-----------------预报精度评定--------------------%
Y_rand=zeros(1,num);
for i=1:num
Y_rand(i)=normrnd(0,sqrt(sigma));%随机项
end
yuce2=X2*fai;
figure
plot(1:num,yuce2,'b',1:num,Y2,'r');
legend('预测值','实测值');
rate2=length(find((yuce2-Y2)./Y2*100<10))/num*100;%合格率
DC2=1-sum((Y2-yuce2).^2)./sum((Y2-mean(Y2)).^2);%确定性系数
|
|