Exercise - Generation of synthetic hydrological series for the Po River at Piacenza This is a possible sequence of instructions in R. The list of instructions is uncommented. It only serves as a guidance. The list assumes that the data are already uploaded in a vector called: "po". Therefore, please note that the data need to be uploaded accordingly. plot(po,type="l") mean(po) min(po) max(po) po1=log(po) plot(po1,type="l") po2=array(po1,dim=c(365,86)) mu=rep(0,365) sigma=rep(0,365) for(i in 1:365) mu[i]=mean(po2[i,]) plot(mu,type="l") for(i in 1:365) sigma[i]=sd(po2[i,]) plot(sigma,type="l") plot(sigma*10,type="l") lines(mu,col="red") po3=(po1-rep(mu,86))/rep(sigma,86) plot(po3,type="l") acf(po1,lag.max=365) # autocorrelation function useful to show the presence of seasonality acf(po1,lag.max=365*87) acf(po3,lag.max=365*87) # plot(acf(po1),acf(po3), type="l",col="black","red") phi1=acf(po3,lag.max=730)\$acf[2] phi1 phi1=acf(po3,lag.max=730)\$acf[1] phi1 phi1=acf(po3,lag.max=730)\$acf[2] epsilon=rep(0,365*86) for(i in 2:(365*86)) epsilon[i]=po3[i]-phi1*po3[i-1] plot(epsilon,type="l") plot(density(epsilon),type="l") acf(epsilon,lag.max=365) acf(epsilon) pr1=ar(po3,order.max=2) # change 2 moving forward 3,4,5.... names(pr1) acf(pr1\$resid[3:(365*86)]) # if changed above, increase also here mean(epsilon) sd(epsilon) epsilon1=rnorm(365*1000,0,sd(epsilon)) plot(epsilon1,type="l") pos1=rep(0,(365*1000)) for(i in 2:(365*1000)) pos1[i]=phi1*pos1[i-1]+epsilon1[i] plot(pos1[1:3650],type="l") acf(pos1[1:3650],lag.max=365) pos2=pos1*rep(sigma,1000)+rep(mu,1000) acf(pos2[1:3650],lag.max=365) pos3=exp(pos2) mean(po) mean(pos3) sd(po) sd(pos3) max(po) max(pos3) plot(density(po),type="l") lines(density(pos3),col="red",lwd=2) # peak lower due to distribution not gaussian # pr1=acf(pos3)\$acf # \$acf for the autocorrelation coefficients pr1=acf(pos3,lag.max=365)\$acf acf(po,lag.max=365) lines(pr1,col="red",lwd=2) # for droughts is very fine