MATLAB: Wind speed ARMA simulation

armasimulatewindspeed

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

Best Answer

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.
Related Question