MATLAB一维自回归预报模型的建立与预报值的复原

MATLAB一维自回归预报模型的建立与预报值的复原

clc;clear all;close all
%% 读取数据
f    = 'NCEP.sst.mean.nc';
lon  = ncread(f,'lon');
lat  = ncread(f,'lat');
time = ncread(f,'time');
sst  = ncread(f,'sst');
 
latq = 25;lonq = 130;%插值点
 
length_1  = length(time);
[LAT LON] = meshgrid(lat,lon);%网格化
%%  插值
sst(sst<-9e36) = 0;
for i = 1:length_1
    sstq(i) = interp2(LAT,LON,squeeze(sst(:,:,i)),latq,lonq);
end 
%%  画图_时间序列曲线
set(0,'defaultfigurecolor','w')
set(gcf,'units','normalized','position',[0.25 0.25 0.5 0.5]);
figure(1);
plot(time,sstq);
datetick('x','yyyymm','keepticks');
xlabel('year/month');
ylabel('sst');
box off;
title('北纬25度，东经130度位置处海表温度时间序列','fontsize',15,'fontweight','bold');
hold on;
%% 去趋势
sst1 = sstq(:);
dsst = detrend(sst1);
trendline = sst1 - dsst;
 
%% 画图_趋势线
plot(time,trendline);
datetick('x','yyyymm','keepticks');
xlabel('years/months');
ylabel('SST(℃)');
legend('时间序列','趋势线');
 
%% 画图_去趋势
figure(2);
set(0,'defaultfigurecolor','w')
set(gcf,'units','normalized','position',[0.25 0.25 0.5 0.5]);
plot(time,dsst,'r-');
datetick('x','yyyymm','keepticks');
xlabel('years/months');
ylabel('SST(℃)');
title('北纬25度，东经130度位置处去趋势后的海表温度时间序列','fontsize',15,'fontweight','bold');
hold on;
%% 年循环，气候态平均
 
timestr = datestr(time,30);
mm = str2num(timestr(:,5:6));
for i = 1:12
    m_index = find(mm == i);
    climate_sst(i) = mean(dsst(m_index));
    circle_sst(m_index) = climate_sst(i);
    dssta(m_index) = dsst(m_index) - climate_sst(i);
end
%% 画图_气候态月平均
 
plot(time,circle_sst,'k--');
datetick('x','yyyymm','keepticks');
xlabel('year/month');
ylabel('SSC(℃)');
legend('去趋势后的海表温度','气候态月平均');
box off;
 
%% 画图_去气候态月平均
 
figure(3);
set(0,'defaultfigurecolor','w')
set(gcf,'units','normalized','position',[0.25 0.25 0.5 0.5]);
plot(time,dssta);
datetick('x','yyyymm','keepticks');
%legend('去趋势、去气候态月平均');
xlabel('year/month');
ylabel('SSC(℃)');
title('北纬25度，东经130度位置处去趋势和气候态月平均后的海表温度时间序列','fontsize',15,'fontweight','bold');
box off;
hold on;
%% 自回归建模
%矩阵xx
n = length_1;
xx = [];
x = [];
for k = 1:5
    for i = k:n-6+k
        x(i-k+1) = dssta(i);
    end
    xx = cat(2,x',xx);
end
%矩阵y
y = dssta(6:n);
y = y';
b = regress(y,xx);
 
for o = 6:length_1
dssta_x(o) = b(1)*dssta(o-1)+b(2)*dssta(o-2)+b(3)*dssta(o-3)+b(4)*dssta(o-4)+b(5)*dssta(o-5);
end
%% 画图_预报值
plot(time,dssta_x);
datetick('x','yyyymm','keepticks');
legend('去趋势、去气候态月平均','预报值');
 
%% 精度检验
eil=dssta_x(6:end)-dssta(6:end);
lamdal=sqrt(sum(eil.^2)/(length(eil)));
lamdat=std(dssta(6:end),1);
g=lamdal/lamdat;