Neural Network (NARX) time series prediction

I am hoping to predict daily soil respiration (Rs) values using 8 input variables by utilizing the Neural Network (NARX). My data are daily averages for 4 years (2014-2018), with NaN values for the observed Rs during winter since Rs measurements are only conducted during growing season. The 8 input variables are complete and have no NaN values. I generated my code below using the neural network tool (nnstart). After running the ‘time-series response’ plot, I notice that NaN values of Rs are also predicted and gap-filled in the plot, however the NaN values in the beginning and end are not predicted (image attached). Is there a mistake I am making here and how can I extract the complete predicted time-series of Rs?

I also tried plotting the final prediction of model (y), but it only contains prediction where observed Rs is not NaN. (image attached). I apprecite your help.

Here is my code:

%% Neural Network

% Step 1: Prepare your data
NN_Variables  = synchronize(TT_Ts_Daily_TR, TT_SM_Daily_TR, TT_Ta_Daily_TR, TT_Energy_Daily_TR, TT_PPT_Daily_TR,Chamber_daily_mean);

input_all  = ([NN_Variables.Ts_Pit,NN_Variables.SM_Pit,NN_Variables.Ta,NN_Variables.PPT,NN_Variables.Rn,NN_Variables.LE,NN_Variables.H,NN_Variables.Gflux]);
Rs_all     = ([NN_Variables.Rs]);

X = tonndata(input_all,false,false);
T = tonndata(Rs_all,false,false);

% Choose a Training Function
% For a list of all training functions type: help nntrain
% 'trainlm' is usually fastest.
% 'trainbr' takes longer but may be better for challenging problems.
% 'trainscg' uses less memory. Suitable in low memory situations.
trainFcn = 'trainscg';  % Scaled conjugate gradient backpropagation.

% Create a Nonlinear Autoregressive Network with External Input
inputDelays = 1:2;
feedbackDelays = 1:2;
hiddenLayerSize = 8;
net = narxnet(inputDelays,feedbackDelays,hiddenLayerSize,'open',trainFcn);

% Choose Input and Feedback Pre/Post-Processing Functions
% Settings for feedback input are automatically applied to feedback output
% For a list of all processing functions type: help nnprocess
% Customize input parameters at: net.inputs{i}.processParam
% Customize output parameters at: net.outputs{i}.processParam
net.inputs{1}.processFcns = {'removeconstantrows','mapminmax'};
net.inputs{2}.processFcns = {'removeconstantrows','mapminmax'};

% Prepare the Data for Training and Simulation
% The function PREPARETS prepares timeseries data for a particular network,
% shifting time by the minimum amount to fill input states and layer
% states. Using PREPARETS allows you to keep your original time series data
% unchanged, while easily customizing it for networks with differing
% numbers of delays, with open loop or closed loop feedback modes.
[x,xi,ai,t] = preparets(net,X,{},T);

% Setup Division of Data for Training, Validation, Testing
% For a list of all data division functions type: help nndivision
net.divideFcn = 'dividerand';  % Divide data randomly
net.divideMode = 'time';  % Divide up every sample
net.divideParam.trainRatio = 70/100;
net.divideParam.valRatio = 15/100;
net.divideParam.testRatio = 15/100;

% Choose a Performance Function
% For a list of all performance functions type: help nnperformance
net.performFcn = 'mse';  % Mean Squared Error

% Choose Plot Functions
% For a list of all plot functions type: help nnplot
net.plotFcns = {'plotperform','plottrainstate', 'ploterrhist', ...
    'plotregression', 'plotresponse', 'ploterrcorr', 'plotinerrcorr'};

% 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);

% Recalculate Training, Validation and Test Performance
trainTargets = gmultiply(t,tr.trainMask);
valTargets = gmultiply(t,tr.valMask);
testTargets = gmultiply(t,tr.testMask);
trainPerformance = perform(net,trainTargets,y);
valPerformance = perform(net,valTargets,y);
testPerformance = perform(net,testTargets,y);

% View the Network
view(net)

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

% Multi-step Prediction
% Sometimes it is useful to simulate a network in open-loop form for as
% long as there is known output data, and then switch to closed-loop form
% to perform multistep prediction while providing only the external input.
% Here all but 5 timesteps of the input series and target series are used
% to simulate the network in open-loop form, taking advantage of the higher
% accuracy that providing the target series produces:
numTimesteps = size(x,2);
knownOutputTimesteps = 1:(numTimesteps-5);
predictOutputTimesteps = (numTimesteps-4):numTimesteps;
X1 = X(:,knownOutputTimesteps);
T1 = T(:,knownOutputTimesteps);
[x1,xio,aio] = preparets(net,X1,{},T1);
[y1,xfo,afo] = net(x1,xio,aio);
% Next the the network and its final states will be converted to
% closed-loop form to make five predictions with only the five inputs
% provided.
x2 = X(1,predictOutputTimesteps);
[netc,xic,aic] = closeloop(net,xfo,afo);
[y2,xfc,afc] = netc(x2,xic,aic);
multiStepPerformance = perform(net,T(1,predictOutputTimesteps),y2)
% Alternate predictions can be made for different values of x2, or further
% predictions can be made by continuing simulation with additional external
% inputs and the last closed-loop states xfc and afc.

% 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'];
view(nets)
[xs,xis,ais,ts] = preparets(nets,X,{},T);
ys = nets(xs,xis,ais);
stepAheadPerformance = perform(nets,ts,ys);

% Deployment
% Change the (false) values to (true) to enable the following code blocks.
% See the help for each generation function for more information.
if (false)
    % Generate MATLAB function for neural network for application
    % deployment in MATLAB scripts or with MATLAB Compiler and Builder
    % tools, or simply to examine the calculations your trained neural
    % network performs.
    genFunction(net,'myNeuralNetworkFunction');
    y = myNeuralNetworkFunction(x,xi,ai);
end
if (false)
    % Generate a matrix-only MATLAB function for neural network code
    % generation with MATLAB Coder tools.
    genFunction(net,'myNeuralNetworkFunction','MatrixOnly','yes');
    x1 = cell2mat(x(1,:));
    x2 = cell2mat(x(2,:));
    xi1 = cell2mat(xi(1,:));
    xi2 = cell2mat(xi(2,:));
    y = myNeuralNetworkFunction(x1,x2,xi1,xi2);
end
if (false)
    % Generate a Simulink diagram for simulation or deployment with.
    % Simulink Coder tools.
    gensim(net);
end


% Extract the output predictions of NN
outputPredictions = (cell2mat(y).*0.0216);

start_date  = datetime('2014-01-01');
end_date    = datetime('2018-12-29');
time_vector = start_date:end_date;

plot(outputPredictions)

[Shreeya]

Hi

Could you explain why you are including data points with NaN prediction in your training ?

Answer 1

[Milan Bansal]

Hi,

As per my understanding you are not getting the expected plot for the "time-series response" plot as the NaN values of "Rs" are also predicted by Neural Network.

When using the synchronize function for preparing the data, the NaN values present in the "Rs" are getting replaced by the interpolated values. Exclude the "Rs" variable in the "synchronize" function to avoid the interpolation and hence randomly predicted values in "time-series response" plot.

Refer to the documentation link to know more about "synchronize" function.

https://www.mathworks.com/help/matlab/ref/timeseries.synchronize.html?s_tid=doc_ta

Refer to MATLAB File Exchange link. This contains information and example code for the similar problem

https://www.mathworks.com/matlabcentral/fileexchange/93000-autoregression-multivariate-prediction-with-a-narx-nn?s_tid=answers_rc2-1_p4_Topic

[Farbod Tabaei]

This worked! Thank you very much.