In tito-am/prsim.tools: Librairie complémentaire à PRSIM. Cette librarie permet de simuler, vérifier et créer les figures pour l'analyse des séries stochastiques.

This notebook is intended to present a preliminary analysis of the first batch of simulations done with PRSIM.

CCF

First of all, let's consider the following excerpt from the Penn State Statistics Website:

The basic problem we’re considering is the description and modeling of the relationship between two time series. In the relationship between two time series ( and ), the series may be related to past lags of the x-series. The sample cross correlation function (CCF) is helpful for identifying lags of the x-variable that might be useful predictors of .

Gatineau

######################################################################
#Tests sur la Gatineau                                               #
#Bassins ? tenir en compte: Cabonga,Baskatong,Maniwaki,Paugan,Chelsea#
######################################################################
rm(list=ls())
load('/media/tito/TIIGE/PRSIM/obs_outaouais.Rdata')
bvs<-names(tests)
#names(stoch_sim)<-bvs

#index<-match(c('Cabonga','Baskatong','Maniwaki','Paugan','Chelsea'),bvs)

sim<-list()
filename<-paste('/media/tito/TIIGE/PRSIM/0.9995/sims_final/stoch_sim_10_outaouais_Kappa_1_9995_LD.Rdata')

load(filename)
#choix_graphe<- readline(prompt="Quel système? (1: Gatineau,2:Outaouais Supérieur,) : ")
choix_graphe<-1
if(choix_graphe==1){
  for(i in c('Cabonga','Baskatong','Maniwaki','Paugan','Chelsea')){
    sim[[i]]<-stoch_sim[[i]]$simulation
  }
} else {

  for(i in c('Dozois','Lac Victoria et lac Granet','Rapide-7','Kipawa','Lac des Quinze','Lac Temiscamingue a Angliers')){
    sim[[i]]<-stoch_sim[[i]]$simulation
  }

}


###


# #define plotting colors
 col_sim<-adjustcolor("#fd8d3c",alpha=0.8)
 col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)
 col_obs <- adjustcolor( "black", alpha.f = 0.2)

# 
# ### plot cross-correlation function
par(mfrow=c(length(sim),length(sim)),mar=c(2,2,2,2))#

### run through each station comtination
#meilleure façon d'avoir un seul xlabel y ylabel c'est de passer par ggplot2.
#https://stackoverflow.com/questions/20163877/how-to-have-a-single-xlabel-and-ylabel-for-multiple-plots-on-the-same-page
bvs<-names(sim)
for(j in 1:length(sim)){
  for(i in 1:length(sim)){
    ### ccf of observations
    data_mat <- matrix(unlist(lapply(sim, "[", , "Qobs")),ncol=length(sim))
    ccf_obs <- ccf(data_mat[,i],data_mat[,j],plot=FALSE)
    ### plot ccfs of observations
    bv_bv_ccf<-paste(bvs[i],'-',bvs[j],sep='')
    plot(ccf_obs$lag,ccf_obs$acf,col=col_obs,type="l",ylim=c(0,1),main=bv_bv_ccf)#,xaxt='n',yaxt='n'

    ### simulated ccf
    ### run through each simulation run
    for(r in 1:10){
      data_mat_sim <- matrix(unlist(lapply(sim, "[", , paste("r",r,sep=""))),ncol=length(sim))
      ccf_sim <- ccf(na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,1],na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,2],plot=FALSE)
      ### add one ccf plot per simulation run
      lines(ccf_obs$lag,ccf_sim$acf,col=col_sim)
      grid (NULL,NULL, lty = 1, col = "grey") 
    }
    ### overplot observations again
    lines(ccf_obs$lag,ccf_obs$acf,col="black",lwd=2)
  }
}

Outaouais Superieur

rm(list=ls())
load('/media/tito/TIIGE/PRSIM/obs_outaouais.Rdata')
bvs<-names(tests)
#names(stoch_sim)<-bvs

#index<-match(c('','Baskatong','Maniwaki','Paugan','Chelsea'),bvs)

sim<-list()
filename<-paste('/media/tito/TIIGE/PRSIM/0.9995/sims_final/stoch_sim_10_outaouais_Kappa_1_9995_LD.Rdata')

load(filename)
#choix_graphe<- readline(prompt="Quel système? (1: Gatineau,2:Outaouais Supérieur,) : ")
choix_graphe<-1
if(choix_graphe==1){
  for(i in c('Dozois','Lac Victoria et Lac Granet','Rapide-7','Rapide-2','Lac des Quinze')){
    sim[[i]]<-stoch_sim[[i]]$simulation
  }
} else {

  for(i in c('Dozois','Lac Victoria et lac Granet','Rapide-7','Kipawa','Lac des Quinze','Riviere Kinojevis')){
    sim[[i]]<-stoch_sim[[i]]$simulation
  }

}


###


# #define plotting colors
 col_sim<-adjustcolor("#fd8d3c",alpha=0.8)
 col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)
 col_obs <- adjustcolor( "black", alpha.f = 0.2)

# 
# ### plot cross-correlation function
par(mfrow=c(length(sim),length(sim)),mar=c(2,2,2,2))#

### run through each station comtination
#meilleure façon d'avoir un seul xlabel y ylabel c'est de passer par ggplot2.
#https://stackoverflow.com/questions/20163877/how-to-have-a-single-xlabel-and-ylabel-for-multiple-plots-on-the-same-page
bvs<-names(sim)
for(j in 1:length(sim)){
  for(i in 1:length(sim)){
    ### ccf of observations
    data_mat <- matrix(unlist(lapply(sim, "[", , "Qobs")),ncol=length(sim))
    ccf_obs <- ccf(data_mat[,i],data_mat[,j],plot=FALSE)
    ### plot ccfs of observations
    bv_bv_ccf<-paste(bvs[i],'-',bvs[j],sep='')
    plot(ccf_obs$lag,ccf_obs$acf,col=col_obs,type="l",ylim=c(0,1),main=bv_bv_ccf)#,xaxt='n',yaxt='n'

    ### simulated ccf
    ### run through each simulation run
    for(r in 1:10){
      data_mat_sim <- matrix(unlist(lapply(sim, "[", , paste("r",r,sep=""))),ncol=length(sim))
      ccf_sim <- ccf(na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,1],na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,2],plot=FALSE)
      ### add one ccf plot per simulation run
      lines(ccf_obs$lag,ccf_sim$acf,col=col_sim)
      grid (NULL,NULL, lty = 1, col = "grey") 
    }
    ### overplot observations again
    lines(ccf_obs$lag,ccf_obs$acf,col="black",lwd=2)
  }
}

Riviere Montreal

rm(list=ls())
load('/media/tito/TIIGE/PRSIM/obs_outaouais.Rdata')
bvs<-names(tests)
#names(stoch_sim)<-bvs

#index<-match(c('','Baskatong','Maniwaki','Paugan','Chelsea'),bvs)

sim<-list()
filename<-paste('/media/tito/TIIGE/PRSIM/0.9995/sims_final/stoch_sim_10_outaouais_Kappa_1_9995_LD.Rdata')

load(filename)
#choix_graphe<- readline(prompt="Quel système? (1: Gatineau,2:Outaouais Supérieur,) : ")
choix_graphe<-1
if(choix_graphe==1){
  for(i in c('Mistinikon','Lady Evelyn','Lower Notch et Indian Chute')){
    sim[[i]]<-stoch_sim[[i]]$simulation
  }
} else {

  for(i in c('Mistinikon','Lady Evelyn','Lower Notch et Indian Chute')){
    sim[[i]]<-stoch_sim[[i]]$simulation
  }

}


###


# #define plotting colors
 col_sim<-adjustcolor("#fd8d3c",alpha=0.8)
 col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)
 col_obs <- adjustcolor( "black", alpha.f = 0.2)

# 
# ### plot cross-correlation function
par(mfrow=c(length(sim),length(sim)),mar=c(2,2,2,2))#

### run through each station comtination
#meilleure façon d'avoir un seul xlabel y ylabel c'est de passer par ggplot2.
#https://stackoverflow.com/questions/20163877/how-to-have-a-single-xlabel-and-ylabel-for-multiple-plots-on-the-same-page
bvs<-names(sim)
for(j in 1:length(sim)){
  for(i in 1:length(sim)){
    ### ccf of observations
    data_mat <- matrix(unlist(lapply(sim, "[", , "Qobs")),ncol=length(sim))
    ccf_obs <- ccf(data_mat[,i],data_mat[,j],plot=FALSE)
    ### plot ccfs of observations
    bv_bv_ccf<-paste(bvs[i],'-',bvs[j],sep='')
    plot(ccf_obs$lag,ccf_obs$acf,col=col_obs,type="l",ylim=c(0,1),main=bv_bv_ccf)#,xaxt='n',yaxt='n'

    ### simulated ccf
    ### run through each simulation run
    for(r in 1:10){
      data_mat_sim <- matrix(unlist(lapply(sim, "[", , paste("r",r,sep=""))),ncol=length(sim))
      ccf_sim <- ccf(na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,1],na.omit(cbind(data_mat_sim[,i],data_mat_sim[,j]))[,2],plot=FALSE)
      ### add one ccf plot per simulation run
      lines(ccf_obs$lag,ccf_sim$acf,col=col_sim)
      grid (NULL,NULL, lty = 1, col = "grey") 
    }
    ### overplot observations again
    lines(ccf_obs$lag,ccf_obs$acf,col="black",lwd=2)
  }
}

tito-am/prsim.tools documentation built on March 17, 2021, 9:09 p.m.