Wind speed prediction using ARIMA model

Dear All,

I am trying to predicte the next 2 hours wind speed of 10-min wind speed reading (12-point ahead forecasting). for that i am trying to compare an ANN-NAR model with ARIMA model. for the last one i am getting problems in the predicted wind speed.

in my code i am using a very simple method which has the following steps:

1- Calculated the Autocorrelation & Partial Autocorrelation functions on the row data in order to:

   A-  see if there is a need for data differencing (Identifiy the d value of the ARIMA model)
   B- try to identify the p,q order of the AR and MA filters respectivly.

2- re-calculate the Autocorrelation & Partial Autocorrelation function on the differenced data in order to see if it changes and to identifiy the correct d-value of the ARIMA model.

3- as this Autocorrelation calculation is time consuming it can be shutdown by the if condition.

4- after reading the results of the "correlation- test" an ARIMA model is created, in this mfile i have created random ARIMA models which need to be tested.

5- after building the ARIMA model the standerd steps are carried out:

   A- estimating the model
   B- calculating the residual (optional at the time being don't care of the results)
   C- testing the residual (optional at the time being don't care of the results)
   D- forecasting the wind speed for the next 2 hours (12- points ahead of 10-min reading)
   E- running the simulation (Monte Carlo simulation) (optional)
   F- comparing the forecasted results with the measured results.

MY QUESTION: Why i am getting a non-realistic forecasting results of wind speed when using the forecast function on the selected (estimated) ARIMA model?

my code is as the following:

if true
  %%introduction for the usage of this mfile
% Creating an ARIMA model in matlab to intgrate with the ANN script to
% compare the different prediction results of those two differnt models in
% the report was called:
% 1- STSFMs-ARIMA
% 2- ANNTSFMs_NAR
clear
clc
tic
%%pre-processing of the wind speed data by doing the ac.f and pac.f
% this part will try to describe the available data which will help use
% understand the best fit model to be used in next steps.
% ------------- Loading the data --------------------
load('C:\Users\00037218\Desktop\Wind Speed prediction\E82_Wind_2_2013.mat')
%----------------------- End ------------------------
% for doing the AutoCorrelation Function and the Partial AutoCorrelation
% Function please use the following functions in matlab:
k_max = round(length(Final_test)/3);% maximum number of lags to be correlated
Test = 0;
if Test
  tic
  [ACF, lags, bounds] = autocorr(Final_test,k_max);
  Feedback_delays = lags((ACF > bounds(1,1)) | (ACF < bounds(2,1)));
  q = Feedback_delays;
  [PACF, Plags, Pbounds] = parcorr(Final_test,k_max);
  P = lags((PACF > Pbounds(1,1)) | (PACF < Pbounds(2,1)));
  figure();
  subplot(2,1,1)
  lineHandles = stem(lags,ACF,'filled','r-o');
  set(lineHandles(1),'MarkerSize',4)
  grid('on')
  xlabel('Lag')
  ylabel('Sample Autocorrelation')
  title('Sample Autocorrelation Function')
  hold('on')

    plot(lags,ACF);hold on; plot(lags,repmat(bounds(1,1),1,max(lags)+1));hold on;plot(lags,repmat(bounds(2,1),1,max(lags)+1)) 
    subplot(2,1,2)
    lineHandles = stem(Plags,PACF,'filled','r-o');
    set(lineHandles(1),'MarkerSize',4)
    grid('on')
    xlabel('Lag')
    ylabel('Sample Partial Autocorrelations')
    title('Sample Partial Autocorrelation Function')
    hold('on')
    plot(Plags,PACF);hold on; plot(Plags,repmat(Pbounds(1,1),1,max(Plags)+1));hold on;plot(Plags,repmat(Pbounds(2,1),1,max(Plags)+1)) 
    toc
end
% if we have a linearly decaying sample ACF indicates a nonstationary
% process (time-serie) for that a diffrencing should be done using the diff
% function in matlab as follows:
if Test
    tic
    Final_test_stationary = diff(Final_test);
    [ACF_D, lags_D, bounds_D] = autocorr(Final_test_stationary,k_max);
    Feedback_delays_D = lags_D((ACF_D > bounds_D(1,1)) | (ACF_D < bounds_D(2,1)));
    q_D = Feedback_delays_D;
    [PACF_D, Plags_D, Pbounds_D] = parcorr(Final_test_stationary,k_max);
    P_D = lags((PACF_D > Pbounds_D(1,1)) | (PACF_D < Pbounds_D(2,1)));
    figure();
    subplot(2,1,1)
    lineHandles = stem(lags_D,ACF_D,'filled','r-o');
    set(lineHandles(1),'MarkerSize',4)
    grid('on')
    xlabel('Lag')
    ylabel('Sample Autocorrelation')
    title('Sample Autocorrelation Function')
    hold('on')
    plot(lags_D,ACF_D);hold on; plot(lags_D,repmat(bounds_D(1,1),1,max(lags_D)+1));hold on;plot(lags_D,repmat(bounds_D(2,1),1,max(lags_D)+1)) 
    subplot(2,1,2)
     lineHandles = stem(lags_D,PACF_D,'filled','r-o');
   set(lineHandles(1),'MarkerSize',4)
   grid('on')
   xlabel('Lag')
   ylabel('Sample Partial Autocorrelations')
   title('Sample Partial Autocorrelation Function')
   hold('on')
   plot(Plags_D,PACF_D);hold on; plot(Plags_D,repmat(Pbounds_D(1,1),1,max(Plags_D)+1));hold on;plot(Plags_D,repmat(Pbounds_D(2,1),1,max(Plags_D)+1)) 
    toc
end

   %%creating the ARIMA-mode
% this part of this mfile will show how to creat an STSFMs-ARIMA model and
% how to adjust the different setting of such a model in matlab
if 0 == 1
P = 1;
D = 1;
Q = 1;
STSFMs_ARIMA = arima(P,D,Q); % the creation of an ARIMA model with the following factors: AR-Factor = 2 (first input), Integration_Factor = 2 (second input), MA_Factor = 2 (third input).
% the built ARIMA model will looks like this:
%     ARIMA(2,2,2) Model:
%     --------------------
%     Distribution: Name = 'Gaussian'
%                P: 4
%                D: 2
%                Q: 2
%         Constant: NaN
%               AR: {NaN NaN} at Lags [1 2]
%              SAR: {}
%               MA: {NaN NaN} at Lags [1 2]
%              SMA: {}
%         Variance: NaN
% now to adjust the settings of such a model (seasonallity, lags ,and so
% on) we can use the following
STSFMs_ARIMA.Seasonality = 12;
% the model now will looks like this:
%     ARIMA(2,2,2) Model Seasonally Integrated:
%     ------------------------------------------
%     Distribution: Name = 'Gaussian'
%                P: 16
%                D: 2
%                Q: 2
%         Constant: NaN
%               AR: {NaN NaN} at Lags [1 2]
%              SAR: {}
%               MA: {NaN NaN} at Lags [1 2]
%              SMA: {}
%      Seasonality: 12
%         Variance: NaN
% as it can be seen from this detailed represntation of the ARIMA model the
% following paramerters can be adjusted:
% 1- D (defrencing Factor)
% 2- Constant (Constant term of the selected model depend on the predicted time series)
% 3- AR (Non-seasonal AutoRegression coefficents)
% 4- SAR (Seasonal AutoRegression Coefficents)
% 5- MA (Non-seasonal MovingAverage coefficents)
% 6- SMA (Seasonal MovingAverage coefficents)
% 7- Seasonality 
% 8- Variance.
% all the previouslly stated parameters could be ajusted in the same way
% the seasonality of the model was changed
STSFMs_ARIMA.AR = [];
STSFMs_ARIMA.MA = [];
Lags = num2cell(nan(1,STSFMs_ARIMA.Seasonality));
STSFMs_ARIMA.SAR = Lags;
STSFMs_ARIMA.SMA = Lags;
else
    STSFMs_ARIMA = arima('SMALags',12,'D',1,'SARLags',12,...
        'Seasonality',12);
    D1 = LagOp({1 -1},'Lags',[0,1]);
    D12 = LagOp({1 -1},'Lags',[0,12]);
    D = D1*D12;
    dFinal_test = filter(D,Final_test);
    STSFMs_ARIMA_1 = arima('SMALags', 12, 'D',1,'SARLags',12,'MALags',1,'ARLags',1,...
        'Seasonality',12);
    STSFMs_ARIMA_Minitab = arima('SMALags', 12,'SMA',0.9852, 'D',1,'SARLags',12,'SAR',0.1264,'MALags',1,'MA',0.0500,'ARLags',1,'AR',0.9688,...
        'Seasonality',12,'Constant',-0.0002256,'Variance',2);
    max_order = 3;
    [aic_matrix,bic_matrix] = ARMA_model_order_tuning(Final_test,max_order);
end
   %%tuning (training or estimating) stage of the built ARIMA model
% for tuning our ARIMA model we will need to use the estimate function in
% matlab with a training time series. as following
Est_STSFMs_ARIMA = estimate(STSFMs_ARIMA,Final_test,'Display','full');
Est_STSFMs_ARIMA_1 = estimate(STSFMs_ARIMA_1,Final_test,'Display','full');
Est_STSFMs_ARIMA_Minitab = estimate(STSFMs_ARIMA_Minitab,Final_test,'Display','full');
% ther reuslts of tuning such model will be shown in matlab as following:
%     ARIMA(2,2,2) Model Seasonally Integrated:
%     ------------------------------------------
%     Conditional Probability Distribution: Gaussian
% 
%                                   Standard          t     
%      Parameter       Value          Error       Statistic 
%     -----------   -----------   ------------   -----------
%      Constant    0.000339604   0.000342367        0.99193
%         AR{1}       -1.07354      0.012271       -87.4857
%         AR{2}      -0.163044     0.0125264       -13.0161
%         MA{1}     0.00845896    0.00237307        3.56456
%         MA{2}      -0.980359    0.00239701       -408.992
%      Variance       0.551794    0.00897504         61.481
% for changing the estimation settings use the optimset as following:
options = optimset('fmincon');
options = optimset(options,'Algorithm','sqp','TolCon',1e-7,'MaxFunEvals',...
    1000,'Display','iter','Diagnostics','on');
   %%check goodness of fit
% now we will try to check of the residuals are normally distributed and
% uncorrelated by using the infer function in matlab and finaly plot the
% resutls by the following code:
res = infer(Est_STSFMs_ARIMA,Final_test);
res_1 = infer(Est_STSFMs_ARIMA_1,Final_test);
res_Minitab = infer(Est_STSFMs_ARIMA_Minitab,Final_test);
% for plotting the results
figure();
subplot(2,2,1)
plot(res./sqrt(Est_STSFMs_ARIMA.Variance));grid on
title('Standardized Residuals')
subplot(2,2,2)
qqplot(res);grid on
subplot(2,2,3)
autocorr(res)
subplot(2,2,4)
parcorr(res)
figure();
subplot(2,2,1)
plot(res_1./sqrt(Est_STSFMs_ARIMA_1.Variance));grid on
title('Standardized Residuals')
subplot(2,2,2)
qqplot(res_1);grid on
subplot(2,2,3)
autocorr(res_1)
subplot(2,2,4)
parcorr(res_1)
figure();
subplot(2,2,1)
plot(res_1./sqrt(Est_STSFMs_ARIMA_Minitab.Variance));grid on
title('Standardized Residuals')
subplot(2,2,2)
qqplot(res_Minitab);grid on
subplot(2,2,3)
autocorr(res_Minitab,108)
subplot(2,2,4)
parcorr(res_Minitab,108)
   %%additional tests on the residuals
% the Durbin-Watson Statistic test with the P-value
[P_value,DW]=dwtest(res,Final_test);
% the Ljung-Box-test
[hLBQ,P_valueLBQ] = lbqtest(res,'lags',[12,24,36,48]);
% Engle's ARCH test is used to identify residual heteroscedasticity
[HARCH,p_valueARCH] = archtest(res,'lags',[12,24,36,48]);
% if the results of the first vlaues of the used test method are 1 means
% that evidence is establisted that the test is positive and the tested
% relation exist yet 0 indicate the the test fail and no evidence of the
% tested relation exists.

    %%Forecasting stage using the Est_STSFMs_ARIMA (where the coefficents are now known)
% for foreacasting use the forecast matlab function as following:
N = 12;% forecast horizon
[Yc,YcMSE,U] = forecast (Est_STSFMs_ARIMA,N);% the second input is the forecast horizon for more details refer to the created report by the author and matlab documentation.
[Yc_1,YcMSE_1,U_1] = forecast (Est_STSFMs_ARIMA_1,N);
[Yc_Minitab,YcMSE_Minitab,U_Minitab] = forecast (Est_STSFMs_ARIMA_Minitab,12);
% for showing the forecasted results use the plot function in matlab as
% following:
figure();
subplot(3,1,1);
plot(Target_wind);grid on;hold on;plot(YcMSE,'-.r');xlabel('Prediction Step Index');ylabel('Wind Speed in m/s');title('Est_STSFMs_ARIMA_12 model for wind speed forecasting');legend('Expected Outputs','ARIMA Predictions');hold off
subplot(3,1,2);
plot(Target_wind);grid on;hold on;plot(Yc_1,'-.r');xlabel('Prediction Step Index');ylabel('Wind Speed in m/s');title('Est_STSFMs_ARIMA_12 model for wind speed forecasting (with non-seasonal terms)');legend('Expected Outputs','ARIMA Predictions');hold off
subplot(3,1,3);
plot(Target_wind);grid on;hold on;plot(Yc_Minitab,'-.r');xlabel('Prediction Step Index');ylabel('Wind Speed in m/s');title('Est_STSFMs_ARIMA_12 model for wind speed forecasting (with non-seasonal terms as in Minitab)');legend('Expected Outputs','ARIMA Predictions');hold off
   %%Monte Carlo simulation of ARIMA
% for simulating different paths of the built model the Monte Carlo
% simulation can be used in Matlab as following:
number_of_paths = 10;
number_of_observation = N;%length(Final_test);
[Y,E,C] = simulate(Est_STSFMs_ARIMA,number_of_observation,'NumPaths',number_of_paths);%length(Final_test));
% for showing the different results of such simulation use the following
% commands:
figure();
subplot(3,1,1);
plot(Y);grid on
title('Simulated Response')
subplot(3,1,2);
plot(E);grid on
title('Simulated Innovations')
subplot(3,1,3);
plot(Y(:,1),'-.r');grid on;hold on; plot(Target_wind);hold off
title('Compare between Simulated and measured wind speed reading')

   %%compare model output and measured output
% this part will try to show the usage of the compare function in matlab
% for doing so please follow those steps:
nx = 1:10;% the model order
sys = ssest(iddata(Final_test),nx);%#ok<*NASGU> % this will give the optimum model order number
%nx = 1; % chnage this value to that given by the previous analysis
sys = ssest(iddata(Final_test),nx);
compare(Target_wind,sys,N);%length(Final_test));%N)
compare(Final_test,sys,length(Final_test));
toc
  end

[Peter Mills]

what is "Final_test"?

[Peter Mills]

what is the function "ARMA_model_order_tuning"? Thank you in advance

Answer 1

[Hang Qian]

To forecast an ARIMA model, we want to provide both a fitted model as well as data. The former only carries model coefficients, but not observations.

Replace

forecast (Est_STSFMs_ARIMA,N)

by

forecast (Est_STSFMs_ARIMA,N,’Y0’,Y)

It will resolve the anomaly of unrealistic values.

Hang Qian

[Aubai]

Thx very much actually that is my prbolem. :)

[Mohsen]

what is Y0 and Y ?

[JiashenTeh]

Hi Hang,

Why when I simulate the ARMA model which I have fitted to my wind speed data , the wind speed I get is negative? this doesn't make any sense to me.

Can you please explain this? Thank you.

using the example above, let's say i do this simulate(Est_STSFMs_ARIMA,N)

and I get negative wind speed..

Answer 2

[Mohsen]

sorry what is the target_wind ? here :plot(Target_wind)

[Aubai]

it is the target value of the wind which Need to be predicted

Answer 3

[sharad kolte]

HI CAN U SEND THE DAT SET OF WIND SPEED,

[Aubai]

sorry mate but i can't

Answer 4

[Lourence Ferrera]

Is your data has seasonality ? How can I code the data to show that it has seasonality and forecast this using SARIMA or ARIMA in matlab?

[Aubai]

Hallo,

try to search for ' Parametric Trend Estimation ' in the help and you will find your answer.

yes my data is seasonal, and you can follow the matlab help and documentation for deseasonalization of time series data.

finally I can suggest the using the Henderson trend-estimation filters family in order to de-seasonalize the data.

Best Regards

[Lourence Ferrera]

if we use henderson trend -estimation filter then we don't need the unit root test for our data right ?

[Aubai]

I am still an electrical engineer so I am not the expert on statistical stuff.

here you can find a link for usage of such test. finally I can say that the unit root test can be used to "The task of the test is to determine whether the stochastic component contains a unit root or is stationary" according to the given reference.

What I am trying to say that:

- Trend Analysis for the D-part of your time series

- Seasonal Analysis for the S-part of your time series

- Stochastic Analysis for the Z-part (component) of your time series.

hope this helps

Answer 5

[Alex]

Hi Aubui,

Could you explain this statement "k_max = round(length(Final_test)/3);% maximum number of lags to be correlated"?

I could not understand the "max no. of lags to be correlated".

Thanks

AD

Answer 6

[Ling Ji]

what does the "ssest()” mean? Is it an self-defined function?