Wind speed ARMA simulation

Dear All,

I have some problem in simulating my wind speed using ARMA model. When I use the ARMA, i obtain negative wind speed which doesn't make sense. Below, I explain my thought process and and work flow.

1) Plot the wind speed and examine it. This is the how the wind speed for a year on an hourly basis looks like.

2) Then, by using the box-jenkins model selection method. I obtain the ARMA model for my wind speed.

   a) Box and Jenkins first step. Determine stationarity of my data. Graph below   shows the autocorrelation (ACF) and partial autocorrelation (PACF) of the data. It seems that my data is not stationary enough. Hence differencing is needed.

   b) Differencing by one time. Plot the ACF and PACF. It seems that both ACF and PACF gradually tail off. Hence, it suggests an ARMA model. Also it suggested p and q lags of no more than 4.

   c) To determine the best combination of p and q lags, I use the Bayesian information criterion (BIC) method.

      The BIC method suggested that p = 4 and q = 4 are best used to model my wind speed data.

3) Test the model adequacy from mathematical view point.

   a) Check that the residual is normally distributed. It seems that it is normally distributed.

   b) Further confirm with box and Jenkins test. Residual error within 10%. QQ plot suggest it is mostly linear, hence suggest individuality of the residual (non-correlated). ACF and PACF suggest that 100% of the residuals stay within the bands. Hence, the p and q lags chosen are suitable.

4) From mathematical view point, it seems all good and well. However the problem comes when I try to ensure that the model has physical meaning in real world. It doesn't make any sense as the wind speed has negative values!

Please help me dear ARMA users!

This is my code:

    if true
      % Original data
      Y = reshape(WindSpeed(1:8736,1:2),8736*2,1);
      plot(Y);

      % Check for stationarity of the original data
      figure
      subplot(2,1,1)
      autocorr(Y,50)
      subplot(2,1,2)
      parcorr(Y,50)

      % Check for stationarity after differencing the data once
      D1 = LagOp({1,-1},'Lags',[0,1]);
      D2 = D1;
      dY = filter(D2,Y);
      figure
      subplot(2,1,1)
      autocorr(dY,50)
      subplot(2,1,2)
      parcorr(dY,50)

      % Use BIC method to determine ARMA lags
      LOGL = zeros(4,4); %Initialize
      PQ = zeros(4,4);
      for p = 1:4
          for q = 1:4
              mod = arima(p,1,q);
              [fit,~,logL] = estimate(mod,Y,'print',false);
              LOGL(p,q) = logL;
              PQ(p,q) = p+q;
           end
      end
      LOGL = reshape(LOGL,16,1);
      PQ = reshape(PQ,16,1);
      [~,bic] = aicbic(LOGL,PQ+1,100);
      BEST = reshape(bic,4,4);
      % At this point, select the minimum BEST values, the row index is the p value    and column index is the q value

      % Diagnostic check, determine whether the model fit properly or not.
      p = 4; q = 4;
      Mdl = arima(p,1,q);
      EstMdl = estimate(Mdl,Y);
      [E,~] = infer(EstMdl,Y);
      % Shows that the residual E is normally distributed
      hist(E); 
      % Box and jenkins test
      figure
      subplot(2,2,1); plot(E./sqrt(EstMdl.Variance)); title('Standardized Residuals')
      subplot(2,2,2); qqplot(E);
      subplot(2,2,3); autocorr(E)
      subplot(2,2,4); parcorr(E)

      % Reality test
      Ysim = simulate(EstMdl,8736,'Y0',Y);
      plot(Ysim);
    end

Thank you very much in advance

[Ashim]

Well, I am not sure what may have caused the problem. I would go a bit deeper into it shortly. I am interested in similar stuff and working with it at the moment. Did you try detrending the wind speed data and then trying the ARMA model. The wind speed that you are getting are in fact the predicted wind speed.

I let you know about my results if I could solve your issue

Answer 1

[Brendan Hamm]

You have a mean zero process with Normal errors and no presample response, so you are essentially starting your prediction with just a shock. If you wish to have the samples be a prediction into the future, you need to specify these presample responses using the 'Y0' argument to the simulate command. I assume you have some presample data which was used to fit the parameters of your model.

[JiashenTeh]

Hi Brendan,

I do not want to predict future values, rather i want to simulate the ARMA model. But i get negative values, which doesn't make any sense.

Please help me!

Thanks!

[Brendan Hamm]

There is nothing inherent in an ARIMA model that guarantees non-negative solutions. For this reason this may be in inappropriate model. One potential suggestion is to apply a log-transformation to your original data first and then attempt an arima model. Of course you still need to follow the methodology after the transformation, but when you apply the inverse transformation you can guarantee non-negativity. I cannot guarantee this will work in your case.

Just a couple extra notes (not that these help solve the issue): Checking normality with the histogram you have is not indicative of normality. You have 2 few bars to even discern what is happening with the data. You should perform some sort of formal hypothesis test, i.e. jbtest or lillietest. If you want to visualize normality use normplot or probplot.

When you call aicbic you tell it you have 100 sample points when it appears you have much more.

[JiashenTeh]

Dear Brendan,

Thank you for pointing out!

So what can we do if we realized that the residual is not normal? Does it indicates that the model is not suitable and we should look for other more suitable model? say for example, i notice that the residual is not normally distributed and the fitted model is ARMA(2,2). So what I can do next considering that i am very sure it is a time series data?