Sarima_LM

clc;clear;
data1=xlsread('data1.xlsx');
num=132;
data=data1(1:num,:);
test1=data1(num+1:end,:);
y = log(data);
T = length(y);
Mdl = arima('Constant',0,'D',1,'Seasonality',12,...
              'MALags',1,'SMALags',12);         
EstMdl = estimate(Mdl,y);
res = infer(EstMdl,y);
stres = res/sqrt(EstMdl.Variance);
%-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
yfit=y-res;
yfit2008=exp(yfit);
y2008=exp(y);
res1=y2008-yfit2008;
res1=res1(14:end);
res1=res1';
Mdl1 = estimate(Mdl,y);
yF1  = forecast(Mdl1,12,'Y0',y);
yF1=exp(yF1);
res2=test1-yF1;
res2=res2';
resinput=[res1,res2];
resinput=resinput';
T = tonndata(resinput,false,false);
trainFcn = 'trainlm';  
feedbackDelays = 1:2;
hiddenLayerSize = 10;
net = narnet(feedbackDelays,hiddenLayerSize,'open',trainFcn);
[x,xi,ai,t] = preparets(net,{},{},T);
% Setup Division of Data for Training, Validation, Testing
net.divideParam.trainRatio = 70/100;
net.divideParam.valRatio = 15/100;
net.divideParam.testRatio = 15/100;
% Train the Network
[net,tr] = train(net,x,t,xi,ai);
% Test the Network
y = net(x,xi,ai);
e = gsubtract(t,y);
performance = perform(net,t,y);
% Plots
% Uncomment these lines to enable various plots.
%figure, plotperform(tr)
%figure, plottrainstate(tr)
%figure, ploterrhist(e)
%figure, plotregression(t,y)
%figure, plotresponse(t,y)
%figure, ploterrcorr(e)
%figure, plotinerrcorr(x,e)

% Closed Loop Network
% Use this network to do multi-step prediction.
% The function CLOSELOOP replaces the feedback input with a direct
% connection from the output layer.
netc = closeloop(net);
netc.name = [net.name ' - Closed Loop'];
view(netc)
[xc,xic,aic,tc] = preparets(netc,{},{},T);
yc = netc(xc,xic,aic);
closedLoopPerformance = perform(net,tc,yc)

% Step-Ahead Prediction Network
% For some applications it helps to get the prediction a timestep early.
% The original network returns predicted y(t+1) at the same time it is
% given y(t+1). For some applications such as decision making, it would
% help to have predicted y(t+1) once y(t) is available, but before the
% actual y(t+1) occurs. The network can be made to return its output a
% timestep early by removing one delay so that its minimal tap delay is now
% 0 instead of 1. The new network returns the same outputs as the original
% network, but outputs are shifted left one timestep.
nets = removedelay(net);
nets.name = [net.name ' - Predict One Step Ahead'];
[xs,xis,ais,ts] = preparets(nets,{},{},T);
ys = nets(xs,xis,ais);
stepAheadPerformance = perform(nets,ts,ys)
%-----------------------------------
y_double=cell2mat(y);%
t_double=cell2mat(t);
%-------------------------------------------------------------------------
yfit1=y_double(1:end-12);%
yfit1=yfit1';
yfit2008=yfit2008(14+2:end);%
yfit1=yfit1+yfit2008;
y1=y2008(14+2:end);
num2=length(yfit1);
mse2= sqrt(sum((y1-yfit1).^2)) ./ num2;
mae2 = mean(abs(y1- yfit1));
rmse2= sqrt(mean((yfit1-y1).^2));
R2 = 1 - (sum((yfit1-y1).^2) / sum((y1-mean(y1)).^2));
disp(['2008-2018年BP神经网络拟合发病率均方误差MSE为：',num2str(mse2)])
disp(['2008-2018年BP神经网络拟合发病率平均绝对误差MAE为：',num2str(mae2)])
disp(['2008-2018年BP神经网络拟合发病率均方根误差RMSE为：',num2str(rmse2)])
disp(['2008-2018年BP神经网络拟合发病率决定系数为：',num2str(R2)])
%----------------------------------------------------------------------
yF11=y_double(end-11:end);%
yF11=yF11';
yF1=yF1+yF11;
figure
plot(yF1,'r--','LineWidth',1.5)
hold on
plot(test1,'b','LineWidth',2)
axis tight
title('2019年SARIMA预测丙肝发病率曲线')
legend('预测值','真实值','Location','NorthWest')
hold off
%------------------------------------------------------------------------------
num1=12;
mse1 = sqrt(sum((test1-yF1).^2)) ./ num1;
mae1 = mean(abs(test1- yF1));
rmse1 = sqrt(mean((yF1-test1).^2));
R1 = 1 - (sum((yF1-test1).^2) / sum((test1- mean(test1)).^2));
disp(['2019年预测发病率均方误差MSE为：',num2str(mse1)])
disp(['2019年预测发病率平均绝对误差MAE为：',num2str(mae1)])
disp(['2019年预测发病率均方根误差RMSE为：',num2str(rmse1)])
disp(['2019年预测发病率决定系数为：',num2str(R1)])
%save filename5.mat