Ajout de variations saisonnières aux séries chronologiques de la vitesse du vent

Question

Suite d'un R blog intéressant et très utile pour simuler les séries temporelles d'une zone inconnue en utilisant ses paramètres de Weibull.

Bien que cette méthode donne une estimation raisonnablement bonne de la série temporelle dans son ensemble, elle souffre beaucoup lorsque nous cherchons des changements saisonniers. Pour tenir compte des variations saisonnières, je veux utiliser les vitesses maximales saisonnières du vent et effectuer la synthèse des séries temporelles de sorte que la distribution annuelle reste constante, c.-à-d. paramètres de forme et d'échelle (valeurs annuelles).

Je souhaite utiliser les vitesses de vent maximales saisonnières pour le code ci-dessous en utilisant 12 vitesses de vent maximales différentes, une pour chaque mois. Cela permettra de plus grandes vitesses de vent à certains mois et plus bas dans d'autres et devrait égaliser les séries chronologiques résultantes.

Le code suit comme ceci:

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

Shape=2; ## Input Shape parameter (yearly). 
Scale=8 ##Calculated Scale Parameter (yearly). 

MaxSpeed<-17 (##yearly) 
## $$$ 12 values of these wind speed one for each month to be used. The resultant time series should satisfy shape and scale parameters $$ ### 
nStates<-16 

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 SpeedSeries 

Quelqu'un peut-il suggérer où et quels changements je dois apporter au code?

Answer 1

Je ne sais pas beaucoup sur la génération de séries temporelles de la vitesse du vent, mais peut-être ces lignes directrices peuvent vous aider à améliorer votre lisibilité code/réutilisabilité:

# 1 Vous voulez probablement avoir une fonction qui va générer un vent temps de vitesse série donné un certain nombre d'observations et une vitesse maximale saisonnière du vent. Alors d'abord essayer de définir votre code à l'intérieur d'un bloc comme celui-ci:

wind_time_serie <- function(nobs, max_speed){ 
    #some code here 
} 

# 2 Ce faisant, s'il semble que certaines parties de votre code sont utiles pour générer des séries chronologiques de la vitesse du vent, mais ne sont pas séries chronologiques de vitesse du vent, essayez de les mettre en fonctions (par exemple, la partie que vous calculez norm_wpdWind, la partie que vous calculez MTMcum, ...). Ensuite, la partie de votre code au début lorsque votre variable globale define devrait disparaître et devenir des arguments par défaut dans les fonctions.

# 4 Évitez d'utiliser les commentaires de fin de ligne lorsque votre ligne est déjà longue et supprimez les demi-colonnes de fin.

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

#Would become this: 
    #Assumes that the initial state is 1 (this must be changed when concatenating days) 
    State<-InitialState<-1 

Ensuite, votre code devrait être plus réutilisable/lisible par d'autres personnes. Vous avez un exemple ci-dessous de ces lignes directrices appliquées à la partie rnorm:

norm_distrib<-function(maxSpeed, states = 16, shape = 2, scale = 8){ 

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

    #Determine Weibull Probability Distribution. 
    wpdWind<-dweibull(WindSpeed,shape=shape, scale=scale) 

    #Convert weibull/Gaussian distribution to normal distribution. 
    return(wpdWind/sum(wpdWind)) 

    } 

    #Plot normal distribution with the max speed you want (e.g. 17) 
    plot(norm_distrib(17),type = "b", ylab= "frequency", xlab = "Wind Speed")