Following up from an

R blog which is interesting and quite useful to simulate the time series of an unknown area using its Weibull parameters.

Although this method gives a reasonably good estimate of time series as a whole it suffers a great deal when we look for seasonal changes. Let's see an example for the certain input it would produce the mean monthly wind speeds that look like this:

*Jan 7.492608
*

Feb 7.059587

March 7.261821

Apr 7.192106

May 7.399982

Jun 7.195889

July 7.290898

Aug 7.210269

Sept 7.219063

Oct 7.307073

Nov 7.135451

Dec 7.315633It can be seen that the variation in wind speed is not that much and in reality, the variation will change from one month to another. If I were to prioritise a certain month say July and June over months of November and December such that the Weibull remains unchanged. How would I do it?

Any lead or advice to make these change in the code listed in the link above would be of great help.

Here is the sample R code.

*MeanSpeed<-7.29 ## Mean Yearly Wind Speed at the site.
*

Shape=2; ## Input Shape parameter.

Scale=8 ##Calculated Scale Parameter.

MaxSpeed<-17

nStates<-16

These are the inputs in the blog, the MeanSpeed is the average annual wind speed at a location that has Shape and Scale parameters as provided. The MaxSpeed is the maximum speed possible over the year.

I would like to have Maxspeed for each month say Maxspeed_Jan, Maxspeed_feb ...till Maxspeed_dec. All with different values. This should be able to reflect the seasonality in the Wind Speed variations across the year.

Then Calculate the following in a certain way that would reflect this variation in the output time series.

*nRows<-nStates;
*

nColumns<-nStates;

LCateg<-MaxSpeed/nStates;

WindSpeed=seq(LCateg/2,MaxSpeed-LCateg/2,by=LCateg) ## Fine the velocity vector-centered on the average value of each category.

##Determine Weibull Probability Distribution.

wpdWind<-dweibull(WindSpeed,shape=Shape, scale=Scale); # Freqency distribution.

plot(wpdWind,type = "b", ylab= "frequency", xlab = "Wind Speed") ##Plot weibull probability distribution.

norm_wpdWind<-wpdWind/sum(wpdWind); ## Convert weibull/Gaussian distribution to normal distribution.

## Correlation between states (Matrix G)

g<-function(x){2^(-abs(x))} ## decreasing correlation function between states.

G<-matrix(nrow=nRows,ncol=nColumns)

G <- row(G)-col(G)

G <- g(G)

##--------------------------------------------------------

## iterative process to calculate the matrix P (initial probability)

P0<-diag(norm_wpdWind); ## Initial value of the MATRIX P.

P1<-norm_wpdWind; ## Initial value of the VECTOR p.

## This iterative calculation must be done until a certain error is exceeded

## Now, as something tentative, I set the number of iterations

steps=1000;

P=P0;

p=P1;

for (i in 1:steps){

r<-P%*%G%*%p;

r<-as.vector(r/sum(r)); ## The above result is in matrix form. I change it to vector

p=p+0.5*(P1-r)

P=diag(p)}

## $$ ----Markov Transition Matrix --- $$ ##

N=diag(1/as.vector(p%*%G));## normalization matrix

MTM=N%*%G%*%P ## Markov Transition Matrix

MTMcum<-t(apply(MTM,1,cumsum));## From the MTM generated the accumulated

##-------------------------------------------

## Calculating the series from the MTMcum

##Insert number of data sets.

LSerie<-52560; Wind Speed every 10 minutes for a year.

RandNum1<-runif(LSerie);## Random number to choose between states

State<-InitialState<-1;## assumes that the initial state is 1 (this must be changed when concatenating days)

StatesSeries=InitialState;

## Initallise----

## The next state is selected to the one in which the random number exceeds the accumulated probability value

##The next iterative procedure chooses the next state whose random number is greater than the cumulated probability defined by the MTM

for (i in 2:LSerie) {

## i has to start on 2 !!

State=min(which(RandNum1[i]<=MTMcum[State,]));

## if (is.infinite (State)) {State = 1}; ## when the above condition is not met max -Inf

StatesSeries=c(StatesSeries,State)}

RandNum2<-runif(LSerie); ## Random number to choose between speeds within a state

SpeedSeries=WindSpeed[StatesSeries]-0.5+RandNum2*LCateg;

##where the 0.5 correction is needed since the the WindSpeed vector is centered around the mean value of each category.

print(fitdistr(SpeedSeries, 'weibull')) ##MLE fitting of SpeedSeriesThe obtained result must resemble the input Scale and Shape parameters. And instead of getting uniform wind speed of each month the variation will reflect the input max wind speeds of each month.

Thank you.