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
1 个评论
Ashim
2015-11-8
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
采纳的回答
Brendan Hamm
2015-3-2
0 个投票
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.
3 个评论
JiashenTeh
2015-3-2
编辑：JiashenTeh
2015-3-2
Brendan Hamm
2015-3-2
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
2015-3-3
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Conditional Mean Models 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
