How can predict multi step ahead using Narnet

I want to predict the future prices, I have used only the daily historical prices as input. I can predict only one step ahead using this code:

clear all;
clear all;
load('prices.mat'); 
set_size = 1413;
targetSeries =prices(1:set_size);  
targetSeries = targetSeries';
targetSeries_train = targetSeries(1:set_size * (4/5)     );
targetSeries_test  = targetSeries(set_size   * (4/5): end);
targetSeries_train = num2cell(targetSeries_train);
targetSeries_test = num2cell(targetSeries_test);
feedbackDelays = 1:4;
hiddenLayerSize = 10;
net = narnet(feedbackDelays, hiddenLayerSize);
net.inputs{1}.processFcns = {'removeconstantrows','mapminmax'};
[inputs,    inputStates, layerStates,targets] = preparets(net,{},{},     targetSeries_train);
[inputs_test,inputStates_test,layerStates_test,targets_test] = preparets(net,{},{},targetSeries_test);
net.trainFcn = 'trainrp';  % Levenberg-Marquardt
net.performFcn = 'mse';  % Mean squared error
net.plotFcns = {'plotperform','plottrainstate','plotresponse', ...
'ploterrcorr', 'plotinerrcorr'};
[net,tr] = train(net,inputs,targets,inputStates,layerStates);
outputs     = net(inputs_test,inputStates_test,layerStates_test);
errors      = gsubtract(targets_test,outputs);
performance = perform(net,targets_test,outputs)
view(net)
%  netc = closeloop(net);
%  [xc,xic,aic,tc] = preparets(netc,{},{},targetSeries_test);
%  yc = netc(xc,xic,aic);
%  perfc = perform(net,tc,yc)

% nets = removedelay(net);
% [xs,xis,ais,ts] = preparets(nets,{},{},targetSeries_test);
% ys = nets(xs,xis,ais);
% stepAheadPerformance = perform(net,ts,ys)

how I can predict multistep ahead, for example the prediction of 10 days in advance.

Thanks in advance.

[John D'Errico]

Argh! Learn how to use the code formatting tool. I'll fix it this once. See that your code is unreadable as it WAS, but after one click of the mouse it is now readable.

[coqui]

Dear friend,

I found the feedbackdelays=1, I want to verify if my presented code is true.

To predict multistep ahead, for example the prediction of 10 days in advance:

% net = narnet(1:1,10);

% net.trainParam.showWindow = false;

% T = tonndata(Y,false,false);

% [Xs,Xi,Ai,Ts] = preparets(net,{},{},T);

% net = closeloop(train(net,Xs,Ts,Xi,Ai));

% Ypred = nan(11,1);

% Ypred = tonndata(Ypred,false,false);

% Ypred(1:1) = T(end:end);

%[xc,xic,aic,tc] = preparets(net,{},{},Ypred);

% Ypred = fromnndata(net(xc,xic,aic),true,false,false);

Thank you very much.

Answer 1

[Shashank Prasanna]

Take a look at the line where I define Ypred here I specify nans for the number of forecasts I want to make. I specify 4 more for pre-samples. Hope this helps. And ofcourse Y is your data Ypred are the forecasts

net = narnet(1:4,10);
net.trainParam.showWindow = false;

T = tonndata(Y,false,false);
[Xs,Xi,Ai,Ts] = preparets(net,{},{},T);
net = closeloop(train(net,Xs,Ts,Xi,Ai));

Ypred = nan(14,1);
Ypred = tonndata(Ypred,false,false);
Ypred(1:4) = T(end-3:end);

[xc,xic,aic,tc] = preparets(net,{},{},Ypred);
Ypred = fromnndata(net(xc,xic,aic),true,false,false);

[coqui]

Hi Shashank,

Thank you very much for your answer. I have some questions:

1- Is Y the input data(prices (1:1413)) or only the test part of the total data (prices set size*1/5)?

2- Supposing that the data (historical values) ranging from 16 January 2000 to 31 December 2010, the Ypred are the prediction of the next 10 days: 1,2,3,4,5,6,7,8,9 and 10 January 2011 ?

Thanks

[coqui]

What is the meaning of delays ?

As feedbackDelays the x past values that explained every current value from the data?

[Shashank Prasanna]

1. The toolbox will do it for you:

http://www.mathworks.com/help/nnet/ug/divide-data-for-optimal-neural-network-training.html

2. Yes

In my experience, you typically use the returns and not the prices (to make it stationary). Anyway I suspect you did your preliminary analysis.

[coqui]

Thank you friend,

I transformed the Prices in returns :(log (Pt/Pt-1)) as inputs, how can reverse the expected returns in prices?

[Shashank Prasanna]

Do the reverse, for continuously-compounded returns:

Pt(i+1) = Pt(i)*exp[RetSeries(i)]

https://en.wikipedia.org/wiki/Continuously_compounded_nominal_and_real_returns

If you have the Financial Toolbox, you can call these functions:

http://www.mathworks.com/help/finance/tick2ret.html

http://www.mathworks.com/help/finance/ret2tick.html

[coqui]

Hi Shashank,

I really appreciate your help.

Thank you very much.

[coqui]

How i can plot the predicted values after training and validation.

I need the code to plot this:

legend ('original targets', 'network predictions','expected outputs')

thanks

[coqui]

I found this code but for narxnet:

figure;
   plot([cell2mat(targetSeries),nan(1,N);
    nan(1,length(targetSeries)),cell2mat(yPred);
    nan(1,length(targetSeries)),cell2mat(targetSeriesVal)]')
legend('Original Targets','Network Predictions','Expected Outputs')

I need your help to plot this figure for narnet.

Thanks.

[omer triyo]

hi, When i try to predict future values, i got same values by increasing number of forecast. for example, i try to predict future values for 365 days. I got same values after about 150 days. What can i do ? Thanks.

Answer 2

[Greg Heath]

%0. a. I do not have MATLAB on this computer, so some of my code comments may need correcting.

% b. I will not complain if you decide to change your mind and accept my answer

% c. Did you mean to begin with close all instead of two clears?

 clear all;
clear all;
load('prices.mat'); 
set_size = 1413;
targetSeries = prices(1:set_size);  
targetSeries = targetSeries';
targetSeries_train = targetSeries(1:set_size * (4/5)     );
targetSeries_test  = targetSeries(set_size   * (4/5): end);}

%1. setsize*4/5 is not an integer

% 2. Why are you trying to avoid the default validation set used to avoid overtraining an overfit net (i.e., more unknown weights than training equations)?

% 3. WARNING: You will have to override the default net.divideFcn = 'dividerand' and associated net.divide... trn/val/tst ratios 0.7/0.15/0.15. Use net.divideFcn = 'dividetrain' with the default ratios 100/0/0. I don't think any other option allows absense of a val set.

%4. HOWEVER, my recommendation is to use net.divideFcn = 'divideblock' and accept the default ratios 0.7/0.15/0.15 . Otherwise is more trouble than it's worth. You will automatically get all three results at once.

 targetSeries_train = num2cell(targetSeries_train);
targetSeries_test = num2cell(targetSeries_test);

% 5. OK, but check out function tonndata

 feedbackDelays = 1:4;

%6. Use the significant delays indicated by the autocorrelation function. For examples search greg narnet nncorr

 hiddenLayerSize = 10;

% 7. Default value of 10 may not work well. However, it is good for preliminary investigation because Numweights = ( 4+1)*10+(10+1)*1= 61 << Numtrainequations = 1130

 net = narnet(feedbackDelays, hiddenLayerSize);
net.inputs{1}.processFcns = {'removeconstantrows','mapminmax'};

% 8. Unnecessary this is a default

 [inputs, inputStates, layerStates,targets] = preparets(net,{},{}, targetSeries_train);
[inputs_test,inputStates_test,layerStates_test,targets_test] = preparets(net,{},{}, targetSeries_test);
net.trainFcn = 'trainrp';  % Levenberg-Marquardt

%9. INCORRECT. 'trainlm' is L-M;

 net.performFcn = 'mse';  % Mean squared error
net.plotFcns = {'plotperform','plottrainstate','plotresponse', 'ploterrcorr', 'plotinerrcorr'};

% 10. Last 3 statements are defaults. Can omit to simplify code

[net,tr] = train(net,inputs,targets,inputStates,layerStates);

% 11. INCORRECT . You did not change the defaults net.divideFcn = 'dividerand' and net.divideRatio = 0.7/0.15/0.15

 outputs     = net(inputs_test,inputStates_test,layerStates_test);
errors      = gsubtract(targets_test,outputs);
performance = perform(net,targets_test,outputs)
view(net)

netc = closeloop(net);

[xc,xic,aic,tc] = preparets(netc,{},{},targetSeries_test);

%%yc = netc(xc,xic,aic);

perfc = perform(net,tc,yc)

%12. Often closing the loop so degrades performance that you have to train netc

 %%nets = removedelay(net);
% %[xs,xis,ais,ts] = preparets(nets,{},{},targetSeries_test);
% %ys = nets(xs,xis,ais);
% %stepAheadPerformance = perform(net,ts,ys)

how I can predict multistep ahead, for example the prediction of 10 days in advance.

% 13. I DON'T UNDERSTAND the use of removedelay. The significant delays of the autocorrelation function indicate how far you can predict ahead with probabilistic certainty. So 1<= FD <= dmax and the values may be nonconsecutive. I think removedelay removes the lower values of FD with, I assume, a decrease in performance. When my computer is available I will investigate.

%Hope this helps.

%Thank you for formally accepting my answer

Greg

[coqui]

Thank you for your answer. I'll try with your code.

[coqui]

How i can plot the predicted values after training and validation.

I need the code to plot this:

legend ('original targets', 'network predictions','expected outputs')

thanks

[Greg Heath]

I don't understand your problem because you have not posted your attempt.

[coqui]

I found this code but for narxnet:

figure;

 plot([cell2mat(targetSeries),nan(1,N);

    nan(1,length(targetSeries)),cell2mat(yPred);
    nan(1,length(targetSeries)),cell2mat(targetSeriesVal)]')
legend('Original Targets','Network Predictions','Expected Outputs')

how we can plot the same figure for narnet?

Thanks

[Greg Heath]

tr=tr

will show the trn/val/tst indices. Therefore you can use the command

hold on

to overlay the separate trn/val/tst results ('b','g','r') over the original data ('k--').

[coqui]

Thank you very much, it's clear.

Please i need your help to answer this question:

http://www.mathworks.com/matlabcentral/answers/142478-how-can-find-the-optimal-delays-and-number-of-hidden-nodes-in-narnet-for-forecasting-task

Thank you in advance.

[Greg Heath]

Search the NEWSGROUP and ANSWERS using combinations of

 greg
 timedelaynet, narnet or narxnet
 nncorr

 Hub, Hmin:dH:Hmax

Hope this helps.

Greg

[MAT-Magic]

Thanks Greg, I am reading the posts in your google group :)

Answer 3

[coqui]

Thank you Greg,

I find these expressions to compute Hub:

1) Hub = -1 + ceil( (Ntrneq-O) / (MXFD*O + O +1) )

2) Hub = -1 + ceil( ( Ntrneq-O)/ (O +1))

The right expression is 1 or 2?

[Greg Heath]

2 is good for fitnet and patternnet

1 is good for narnet(FD)

For timedelaynet replace MXFD with MXID: timedelaynet(ID)

For narxnet use both

[coqui]

Thank you Greg.

As the code right to find optimal delay?

 N = length(y)  %  TargetSeries

 for k = 1:100
    n                = randn(1,N);
    autocorrn        = nncorr(n,n,N-1,'biased');
    sortabsautocorrn = sort(abs(autocorrn));
    M                = floor(0.95*N) ;
    thresh95(k)      = sortabsautocorrn(M);
 end
sigthresh95 = mean(thresh95)
zy = zscore(y,1);
autocorry = ifft( abs(fft(zy)).^2 )/N ;

 lags        = 0:N;

 siglags     = find(abs(autocorry(N:end))>sigthresh95);

THANKS

[Greg Heath]

Looks OK

[coqui]

Thank you Greg, I count on your high knowledge in the field.

When we can use this code:(N-1 instead of N)

    M                = floor(0.95*(2*N-1))  

    lags        = 0:N-1;

and

 siglags     = -1+find(abs(autocorrzt(N:end))>sigthresh95);

instead of

 siglags     = find(abs(autocorry(N:end))>sigthresh95);

Thanks

[coqui]

By using the above program (to obtain optimal feedback delays),I found several outputs, where is the output which represents the optimal lag.

I need an urgent response please.

Thanks in advance

[Greg Heath]

That should give you all of the lags associated with absolute values of correlations above the 95% confidence level for noise.

Try to use the smallest subset of the smallest lags that will yield your training goal. For example, if the significant lags are 1 3 4 7 ..., I would first try [ 1 3 ]. If unsuccessful, add 4 and 7 etc.

P.S. I just posted a NEWSGROUP example of finding the significant autocorrelation lags for the simplenar_dataset.

[coqui]

Thanks Greg, I have found siglags=1 , it means that the optimal lag equals to one???

minmax(autocorrY) = [ -0.5140 1 ]

[coqui]

I have used these code, the last line code is right?

   N = length(y)  %  TargetSeries

 for k = 1:100
    n                = randn(1,N);
    autocorrn        = nncorr(n,n,N-1,'biased');
    sortabsautocorrn = sort(abs(autocorrn));
    M                = floor(0.95*N) ;
    thresh95(k)      = sortabsautocorrn(M);
 end

 sigthresh95 = mean(thresh95)

 zy = zscore(y,1);

 autocorry = ifft( abs(fft(zy)).^2 )/N ;

 lags        = 0:N;

 siglags     = find(abs(autocorry(N:end))>sigthresh95);

[Greg Heath]

This is incorrect. The length of the output of NNCORR is 2*N-1 whereas the length of the fft method is only N.

BOTTOM LINE :

Use the same method to calculate sigthresh95 and siglags.

[coqui]

I have remplaced N-1 by 2*N-1 but I meet this error message: Error using nncorr (line 69) MAXLAG should be less than the number of timesteps of A and B.

 N = length(y)  %  TargetSeries (3000 values)
 for k = 1:100
    n                = randn(1,N);
    autocorrn        = nncorr(n,n,(2*N-1),'biased');
    sortabsautocorrn = sort(abs(autocorrn));
    M                = floor(0.95*N) ;
    thresh95(k)      = sortabsautocorrn(M);
 end
 sigthresh95 = mean(thresh95)
 zy = zscore(y,1);
 autocorry = ifft( abs(fft(zy)).^2 )/N ;
 lags        = 0:N;
 siglags     = find(abs(autocorry(N:end))>sigthresh95);

[Greg Heath]

No.

Your previous formula for autocorrn was correct. whereas your previous formula for M was incorrect

[coqui]

Thanks Greg I'll try again.

[coqui]

I have found siglags=1, what does it means?

Thanks

[Greg Heath]

I don't believe there is just one significant lag. Please post your data in *.m or *.txt.

[coqui]

Thank you very much greg this is my datasample:

https://www.dropbox.com/s/gik6nfuduvbukz6/DATA.txt?oref=e

[Greg Heath]

 There are 1752 significant lags. The first 20 are given below

 [ O N ]   = size(t)        % [ 1 3263 ]
 zt = zscore(t',1)';
 L = floor(0.95*(2*N-1));
for i=1:100
    zn               = zscore(randn(1,N),1);
    autocorrn        = nncorr(zn,zn,N-1,'biased');
    sortabsautocorrn = sort(abs(autocorrn));
    thresh95(i)      = sortabsautocorrn(L);
end
 sigthresh    = mean(thresh95) % 0.025642
 stdsigthresh = std(thresh95)  % 0.00063742

 autocorrt = nncorr(zt,zt,N-1,'unbiased');
 siglags   = -1+find(abs(autocorrt(N+1:2*N-1))>sigthresh);
 Nsig      = size(siglags)  % 1752
 FD        = siglags(2:21)  
 %    1    2    3    4    5   6    7    8   9  11 
 %   12  13  14  15  16  20  21  25  26 31

 Ideally, choose the smallest subset of the smaller values that will yield an acceptable result.

Hope this helps.

Greg

[coqui]

Thank you very much Greg.

Yes, I have found the following first 28 significant lags:

 %    1    2    3    4    5   6    7    8   9  11 
 %   12  13  14  15  16  20  21  25  26  31  32  35  36  37  38  39  40

I'll test with these different lags the most accurate prediction results.

Thanks

[coqui]

The data represents all data under study (train sample+validationsample+testsample).

to identify optimal lag, we must apply the code on all data or only training sample? Thanks

[Greg Heath]

With a data set of your size, determine parameters with the training set and validate with the validation set. Verify consistency using multiple designs with different initial weights. If you wish, you can use the multiplicity to estimate standard deviations on error estimates.

For stability I try to minimize the number of delays, size of the delays and number of hidden nodes.

When all is said and done, obtain an UNBIASED estimate using the test set which, in no way, is used to determine the design.

Hope this helps.

Greg

[coqui]

thanks greg, in this case we must use only training set or (training+validation)set to define optimal lag?

[Greg Heath]

The basic reason for the trn/val/tst division is to make sure that the final design and it's estimates are as unbiased w.r.t. the details of the individual design (training+val) data as possible.

Ideally, you would like Ntrn to be large enough to make reliable estimates of the correlation functions. The most practical approach is to use crossvalidation or bootstrapping estimation to estimate reliability.

Hope this helps.

Greg

[coqui]

Thank you Greg.

Yes, I have found the following first 28 significant lags:

 %    1    2    3    4    5   6    7    8   9  11 
 %   12  13  14  15  16  20  21  25  26  31  32  35  36  37  38  39  40

by employed the design sample(training+val) but I want to know how we can define the optimal one from the 28 possiblities.

Based on mse of test???

[Greg Heath]

My advice:

1. Use the training subset to choose the significant lags. Keep validation predictions as unbiased as possible.

2. There is a formal way: Use the method of partial correlations that chooses lags one-by-one taking into account previously selected correlations. However, this method will yield gaps. I have never tried this.

3. You could choose a subset of the significant correlations by ordering the magnitude. However, that will usually also leave gaps; Similarly, ordering significant correlation lags also tends to have gaps (e.g., 9,11 and 17,19).

4. I prefer to have

 a. No gaps
 b. As small a continuous subset as possible containing, possibly, non-significant correlations.
 c. Trial and error might find me first trying 1:4 and, if unsuccessful, maybe trying 1:8 (even if the correlation at one or more of the intermediate lags is not significant).

Your choice.

Hope this helps.

Greg

[coqui]

Thank you Greg.

As there are any formal way. I have chosen the first: 1. Use the training subset to choose the significant lags.

The design definition is linked to training sample results and no to the test sample.

So, the optimal lags found in my case is 31 lags.

what do you think Greg?

[Greg Heath]

Start with a subset containing the first few. If that doesn't work, add more.

If you plot the autocorreltion function, the two threshold lines, and color the points outside of the lines red, you will get a much better idea of what is going on. Similarly for crosscorrelations.

[EanX]

To obtain 10 steps-ahead forecast, I have to use an empty cell array as input to closed loop net or a cell array of NaN? I suppose that both methods will work so a tried this code:

clearvars; clc; close all;
% test with a simple dataset
Y=simplenar_dataset;

net = narnet(1:4,3);
net.trainParam.showWindow = false;
T = Y(1:80);%tonndata(Y,false,false);

% use only first 80 samples to train and the remaining 20
% to test predictions
Tpred=Y(81:end);

[Xs,Xi,Ai,Ts] = preparets(net,{},{},T);

% is necessary to add Xi and Ai to closeloop to have xic1 and aci1
[net, xic1, aic1] = closeloop(train(net,Xs,Ts,Xi,Ai),Xi,Ai);

% NaN method, require to set first 4 (equal to delay) values 
Ypred = nan(14,1);
Ypred = tonndata(Ypred,false,false);
Ypred(1:4) = T(end-3:end);

% empty cells method
Xc1=cell(1,10);

[xc,xic,aic,tc] = preparets(net,{},{},Ypred);

% prediction resulting from the two methods does not agree
% empty cell does not work but perhaps I'm missing something?

% empty cells method
Ypred1 = fromnndata(net(Xc1,xic1,aic1),true,false,false);

Ypred = fromnndata(net(xc,xic,aic),true,false,false); % NaN method

% plot results
plot([Ypred Ypred1 cell2mat(Tpred(1:10))']);
legend('Pred-NaN','Pred-EmptyCells','Target');

But I obtain different results. What am I missing? Thanks.